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Abstract 



We introduce the Reduced-Rank Hidden Markov Model (RR-HMM), a generalization of HMMs that 
can model smooth state evolution as in Linear Dynamical Systems (LDSs) as well as non-log-concave 
predictive distributions as in continuous-observation HMMs. RR-HMMs assume an m-dimensional latent 
state and n discrete observations, with a transition matrix of rank k < m. This implies the dynamics 
evolve in a fc-dimensional subspace, while the shape of the set of predictive distributions is determined 
by to. Latent state belief is represented with a fc-dimensional state vector and inference is carried out 
entirely in K fe , making RR-HMMs as computationally efficient as fc-state HMMs yet more expressive. 
To learn RR-HMMs, we relax the assumptions of a recently proposed spectral learning algorithm for 
HMMs [1| and apply it to learn fc-dimensional observable representations of rank-fc RR-HMMs. The 
algorithm is consistent and free of local optima, and we extend its performance guarantees to cover the 
RR-HMM case. We show how this algorithm can be used in conjunction with a kernel density estimator to 
efficiently model high-dimensional multivariate continuous data. We also relax the assumption that single 
observations are sufficient to disambiguate state, and extend the algorithm accordingly. Experiments on 
synthetic data and a toy video, as well as on a difficult robot vision modeling problem, yield accurate 
models that compare favorably with standard alternatives in simulation quality and prediction capability. 

1 Introduction 

Models of stochastic discrete-time dynamical systems have important applications in a wide range of fields. Hidden Markov 
Models (HMMs) |2| and Gaussian Linear Dynamical Systems (LDSs) [3] are two examples of latent variable models 
of dynamical systems, which assume that sequential data points are noisy, incomplete observations of a latent state that 
evolves over time. HMMs model this latent state as a discrete variable, and represent belief as a discrete distribution over 
states. LDSs on the other hand model the latent state as a set of real-valued variables, are restricted to linear transition and 
observation functions, and employ a Gaussian belief distribution. The distributional assumptions of HMMs and LDSs also 
result in important differences in the evolution of their belief over time. The discrete state of HMMs is good for modeling 
systems with mutually exclusive states that can have completely different observation signatures. The joint predictive 
distribution over observations is allowed to be non-log-concave when predicting or simulating the future, leading to what 
we call competitive inhibition between states (see Figure [3] below for an example). Competitive inhibition denotes the 
ability of a model's predictive distribution to place probability mass on observations while disallowing mixtures of those 
observations. Conversely, the Gaussian joint predictive distribution over observations in LDSs is log-concave, and thus does 
not exhibit competitive inhibition. However, LDSs naturally model smooth state evolution, which HMMs are particularly 
bad at. The dichotomy between the two models hinders our ability to compactly model systems that exhibit both competitive 
inhibition and smooth state evolution. 

We present the Reduced-Rank Hidden Markov Model (RR-HMM), a smoothly evolving dynamical model with the ability to 
represent nonconvex predictive distributions by relating discrete-state and continuous-state models. HMMs can approximate 
smooth state evolution by tiling the state space with a very large number of low-observation-variance discrete states with a 
specific transition structure. However, inference and learning in such a model is highly inefficient due to the large number 
of parameters, and due to the fact that existing HMM learning algorithms, such as Expectation Maximization (EM) Q, 
are prone to local minima. RR-HMMs allow us to reap many of the benefits of large-state-space HMMs without incurring 
the associated inefficiency during inference and learning. Indeed, we show that all inference operations in the RR-HMM 
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Figure 1: (A) The graphical model representation of an RR-HMM. l t denotes the fc-dimensional state vector, ht the Tri- 
dimensional discrete state, and x t the discrete observation. The distributions over h t and l t+ i are deterministic functions of 
It- (B) An illustration of different RR-HMM parameters and the spaces and random variables they act on. (C) Projection of 
sets of predictive distributions of a rank 3 RR-HMM with 10 states, and a 3-state full-rank HMM with similar parameters. 



can be carried out in the low-dimensional space where the dynamics evolve, decoupling their computational cost from the 
number of hidden states. This makes rank-k RR-HMMs (with any number of states) as computationally efficient as te- 
state HMMs, but much more expressive. Though the RR-HMM is in itself novel, its low-dimensional M. k representation is 
related to existing models such as Predictive State Representations (PSRs) [4|, Observable Operator Models (OOMs) 151 , 
generalized HMMs |6|, and weighted automata |7 8|, as well as the the representation of LDSs learned using Subspace 
Identification [9 1. These and other related models and algorithms are discussed further in Section|5j 

To learn RR-HMMs from data, we adapt a recently proposed spectral learning algorithm by Hsu, Kakade and Zhang fl] 
(henceforth referred to as HKZ) that learns observable representations of HMMs using matrix decomposition and regression 
on empirically estimated observation probability matrices of past and future observations. An observable representation of 
an HMM allows us to model sequences with a series of operators without knowing the underlying stochastic transition and 
observation matrices. The HKZ algorithm is free of local optima and asymptotically unbiased, with a finite-sample bound 
on L\ error in joint probability estimates from the resulting model. However, the original algorithm and its bounds assume 
(1) that the transition model is full -rank and (2) that single observations are informative about the entire latent state, i.e. 
1-step observability. We show how to generalize the HKZ bounds to the low-rank transition matrix case and derive tighter 
bounds that depend on fc instead of m, allowing us to learn rank- A; RR-HMMs of arbitrarily large m in 0(Nk 2 ) time, where 
N is the number of samples. We also describe and test a method for circumventing the 1-step observability condition by 
combining observations to make them more informative. A version of this learning algorithm can learn general PSRs ifTOl 
though our error bounds don't yet generalize to this case. 

Experiments show that our learning algorithm can recover the underlying RR-HMM in a variety of synthetic domains. 
We also demonstrate that RR-HMMs are able to compactly model smooth evolution and competitive inhibition in a clock 
pendulum video, as well as in real-world mobile robot vision data captured in an office building. Robot vision data (and, 
in fact, most real-world multivariate time series data) exhibits smoothly evolving dynamics requiring multimodal predictive 
beliefs, for which RR-HMMs are particularly suited. We compare performance of RR-HMMs to LDSs and HMMs on 
simulation and prediction tasks. Proofs and details regarding examples are in the Appendix. 



2 Reduced-rank Hidden Markov Models 

Let ht £ 1 , . . . , m denote the discrete hidden states of an HMM at time t, and x t € 1 , . . . , n denote the discrete observations. 
Let T e R mxm be the state transition probability matrix with ly = Pr[h t +i — i \ h t = j]. O e R nxm is the observation 
probability matrix such that = Pr[x t = i \ h t — j], and 7? S 1™ is the initial state distribution with 7?j = Pr[/ii = i]. 

Let h t € M. m denote the system's belief, i.e. a distribution over hidden states at time t. If we use a to denote the i th 
column of the identity matrix, then h t is equivalent to the conditional expectation of e/ lt , with the conditioning variables 
clear from context. In addition to the standard HMM notation, assume T has rank k and let T = RS where R € R mxk 
and S £ R kxm . This implies that the dynamics of the system can be expressed in R k rather than M. m . By convention, we 
think of S as projecting the 77i-dimensional state distribution vector to a fc-dimensional state vector, and R as expanding 
this low-dimensional state back to an m-dimensional state distribution vector while propagating it forward in time. One 
possible choice for R and S is to use any fc independent columns of T as the columns of R, and let the columns of S contain 
the coefficients required to reconstruct T from R, tho ugh other choices are possible (e.g. using SVD). Also assume for now 



that m < n (we relax this assumption in Section 3.4 1. We denote the fc-dimensional projection of the hidden state vector h 
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as It, which is simply a vector of real numbers rather than a stochastic vector. We assume the initial state distribution lies 
in the low dimensional space as well, i.e. tt — Rifi for some vector 7?; G R k . Figure [TJa) illustrates the graphical model 
corresponding to an RR-HMM. Figure [TjB) illustrates some of the different RR-HMM parameters and the spaces they act 
on. 

To see how the probability of a sequence can be computed using these parameters, define 

A X = RS diag(O ai i, ...,0 Xjm ) 

W X = S dmg{0 X!l) ...,0 X!m )R 
Tlx W x = SR. With these definitions, the joint probability of xi, Xt, can written 

Pr[xi,...,x t ] = \ x m A Xt ...A Xl n (2.1a) 

= l T m RW Xt ...W Xl n (2.1b) 

The latter parametrization casts a rank-fc RR-HMM as a fc-dimensional PSR or transformed PSR ifTTI . Inference can be 
carried out in 0(Nk 2 ) time in this representation. However, since every HMM is trivially a PSR, this leads to the question 
of how expressive rank-fc RR-HMMs are in comparison to fc-state full-rank HMMs. The following example is instructive. 

2. 1 Expressivity of RR-HMMs 

We describe a rank-fc RR-HMM whose set of possible predictive distributions is easy to visualize and describe. Consider 
the following rank 3 RR-HMM with 10 states and 4 observations. The observation probabilities in each state are of the form 

Oi,. = [piqt pi(l~qi) (l-pi)qi (1 ~Pi)(l ~ Qi)] 

for some < Pi,qi < 1, which can be interpreted as 4 discrete observations, factored as two binary components which are 
independent given the state. T and pi , qi are chosen to place the vertices of the set of possible predictive distributions on 
evenly spaced points along a circle in (p, q) -space: 

Tij = (l/2m) [2 + sin (2iri/m) sin {2irj /m) + cos (27ri/m) cos (27rj'/m)] 
Pi = [sin(27ri/m) + 1] /2 
qi — [cos(27ti/to) + 1] /2 

We plot the marginal probability of each component of the observation, ranging across all achievable values of the latent 
state vector for the m = 10 case (Figure[TJC)), yielding a 10-sided polygon as the projection of the set of possible predictive 
distributions. These distributions are the columns of T t O. We also plot the corresponding marginals for the m — 3 full- 
rank HMM case to yield a triangular set. More generally, from a fc-state HMM, we can get at most a fc-sided polygon for 
the set of possible predictive distributions. 

The above example illustrates that rank-fc RR-HMMs with m states can model sets of predictive distributions which full- 
rank HMMs with less than m states cannot express. However, as we shall see, inference in rank-fc RR-HMMs of arbitrary 
to is as efficient as inference in fc-state full -rank HMMs. This implies that the additional degrees of freedom in the RR- 
HMM's low-dimensional parameters and state vectors buy it considerable expressive power. Since RR-HMMs are also 
related to PSRs as pointed out in the previous section, and since our learning algorithm will be shown to be consistent for 
estimating PSRs (though we have finite-sample guarantees only for the RR-HMM case), it is also instructive to examine the 
expressivity of PSRs in general. We refer the reader to Jaeger (2000) [5 1 and James et. al. (2004) 1 12 1 for more on this. 

3 Learning Reduced-Rank HMMs 

In a full-rank HMM, the maximum likelihood solution for the parameters {T, O} can be found through iterative techniques 
such as expectation maximization (EM) fl3l . EM, however, is prone to local optima and does not address the model 
selection problem. HMM model selection algorithms that avoid local minima (e.g. [ 14|) are better but still not guaranteed 
to return anything close to optimal as data increases, and is slow beyond a certain state space magnitude. Moreover, in 
learning RR-HMMs we face the additional challenge of learning the factors of its low-rank transition matrix. We could use 
EM to estimate T followed by (or combined with) matrix factorization algorithms such as Singular Value Decomposition 
(SVD) ifTSIl or Non-negative Matrix Factorization (NMF) lfl6l . This approach has several drawbacks. For example, if the 
noisy estimate of a low-rank transition matrix is not low-rank itself, SVD could cause negative numbers to appear in the 
reconstructed transition matrix. Also, algorithms for NMF are only locally optimal, and NMF is overly restrictive in that it 



so that A x € E mxm , and define 

so that W x € R kxk . Also let W = 
using eitherj^lj.} or {W x } as 
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constrains its factor matrices to be non-negative, which is unnecessary for our application since low-rank transition matrices 
may have negative numbers in their factors R and S. 

An alternative approach, which we adopt, is to learn an asymptotically unbiased observable representation of an RR-HMM 
directly using SVD of a probability matrix relating past and future observations. This idea has roots in subspace identifi- 
cation 12] El and multiplicity automata QE]|6] as well as the PSR/OOM literature ll5l [T8ll and was recently formulated in 
a paper by Hsu, Kakade and Zhang 1 1 1 for full-rank HMMs. We use their algorithm, extending its theoretical guarantees 
for the low-rank HMM case where the rank of the transition matrix T is k < m. Computationally, the only difference in 



our base algorithm (before Section 3.4 1 is that we learn a rank k representation instead of rank m. This allows us learn 
much more compact representations of possibly large-state-space real-world HMMs, and greatly increases the applicability 
of the original algorithm. Even when the underlying HMM is not low-rank, we can examine the singular values to tune 
the complexity of the underlying RR-HMM, thus providing a natural method for model selection. We present the main 
definitions, the algorithm and its performance bounds below. Detailed versions of the supporting proofs and lemmas can be 
found in the Appendix. 



3.1 The Algorithm 

The learning algorithm depends on the following vector and matrix quantities that comprise properties of single observa- 
tions, pairs of observations and triples: 

[Pi]i = Pi[xi = i] 

[P2,l]i,j = Pr[x 2 = i,X!= j] 

[Pj,x,i] = P*[x3 = i,X2 = x,X! = j] for x = 1, . . . , n 

Pi £ R" is a vector, P2.1 £ M" x ™ and Ps :X ,i G R nx " are matrices. These quantities are closely related to matrices 
computed in algorithms for learning OOMs [5], PSRs 1 18 1 and LDSs using subspace identification (Subspace ID) [9|. They 
can be expressed in terms of HMM parameters (for proofs see the Appendix: Lemmas [8] and [9] in Section [7!T) : 

PJ = P m T diag(^)0 T 
P 2 ,i = OT diag(7r)0 T 
Pw = OA x T diag(7r)0 T 

Note that P21 and P3 >x ,i both contain a factor of T and hence are both of rank k for a rank-fc RR-HMM. This property 
will be important for recovering an estimate of the RR-HMM parameters from these matrices. The primary intuition is 
that, when projected onto an appropriate linear subspace, Pz. x ,i is linearly related to P 2 .i through a product of RR-HMM 
parameters. This allows us to devise an algorithm that 

1. estimates P2,i and P^. Xt \ from data, 

2. projects them to an appropriate linear subspace computed using SVD, 

3. uses linear regression to estimate the RR-HMM parameters (up to a similarity transform) from these projections. 

Specifically, the algorithm attempts to learn an observable representation of the RR-HMM using a matrix U £ E nxfc such 
that U T OR is invertible. An observable representation is defined as follows. 

Definition 1 The observable representation is defined to be the parameters b\, 6oo, {B x }™ =1 such that: 

bi = U J Pi (3.1a) 
boo = (PliU) + Pi (3.1b) 
B x = (U T P 3 , X , 1 )(U T P 2 , 1 )+ forx = l,...,n (3.1c) 



For the RR-HMM, note that the dimensionality of the parameters is determined by k, not m: 61 S M. k , b^ £ R k and 
Vx B x £ R kxk . Though these definitions seem arbitrary at first sight, the observable representation of the RR-HMM 
is closely related to the true parameters of the RR-HMM in the following manner (see Lemma [9] in the Appendix for the 
proof): 

1. b x = (U T OR)iri = (U t O)it, 

2. bl = l r m R(U r OR)-\ 
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3. For all a; = l,...,n : B x = {U T OR)W x {U T OR)~ 1 



Hence B x is a similarity transform of the RR-HMM parameter matrix W x — S &ia,g(0 X: .)R (which, as we saw earlier, 

allows us to perform RR-HMM inference), and 61 and are the corresponding linear transformations of the RR-HMM 
initial state distribution and the RR-HMM normalization vector. Note that (U T OR) must be invertible for these relationships 
to hold. Together, the parameters b^b^ and B x for all x comprise the observable representation of the RR-HMM. Our 
learning algorithm will estimate these parameters from data. The algorithm for estimating rank-fc RR-HMMs is equivalent 
to the spectral HMM learning algorithm of HKZ [1] for learning fc-state HMMs. Our relaxation of their conditions (e.g. 
HKZ assume a full-rank transition matrix, without which their bounds are vacuous), and our performance guarantees for 
learning rank-fc RR-HMMs, show that the algorithm leams a much larger class of fc-dimensional models than the class of 
fc-state HMMs. 

LEARN-RR-HMM(fc, N) The learning algorithm takes as input the desired rank fc of the underlying RR-HMM rather 
than the number of states m. Alternatively, given a singular value threshold the algorithm can choose the rank of the HMM 
by examining the singular values of P 2 ,i m Step 2. It assumes that we are given N independently sampled observation 
triples (x%, x%, X3) from the HMM. In practice, we can use a single long sequence of observations as long as we discount 
the bound on the number of samples based on the mixing rate of the HMM (i.e. (1 — the second eigenvalue of T)), in which 
case 7r must correspond to the stationary distribution of the HMM to allow estimation of P\. The algorithm results in an 
estimated observable representation of the RR-HMM, with parameters b\, b^, an d B x for x — 1, . . . , n. The steps are 
briefly summarized here for reference: 

1. Compute empirical estimates Pi, P2,i, Ps,x,\ of Pit P%1, P$,x,i (f° r x = 1, n). 

2. Use SVD on Pi,\ to compute U, the matrix of left singular vectors corresponding to the fc largest singular values. 

3. Compute model parameter estimates: 

(a) $1 = C/ T Pi, 

(b) = (P 2 T 1 C/)+P 1 , 

(c) B x = U T P 3 , X , 1 (U T P 2A )+ (forx= l,...,n) 

We now examine how we can perform inference in the RR-HMM using the observable representation. For this, we will 
need to define the internal state b t - Just as the parameter b\ is a linear transform of the initial RR-HMM belief state, b t is a 



linear transform of the belief state of the RR-HMM at time t (Lemma [T0| in Section 7.1 of the Appendix): 

b t = {U r OR%{x x .. t - X ) = (U 1 '0)h t (x 1:t -i) 

This internal state b t can be updated to condition on observations and evolve over time, just as we can update l t for RR- 
HMMs and h t for regular HMMs. 

3.2 Inference in the Observable Representation 

Given a set of observable parameters, we can predict the probability of a sequence, update the internal state b t to perform 



filtering and predict conditional probabilities as follows (see Lemma 10 in the Appendix for proof): 



Predict sequence probability: Pr[xi, . . . , x t ] — ^bJ X! B Xt . . . B Xl b\ 
Internal state update: bt+i — - bt - 

Conditional probability of x t given Xi-t-i: Pr[ir t | Xi : t-x] 



T, I b T ao B J: b t 



Estimated parameters can, in theory, lead to negative probability estimates. These are most harmful when they cause 
the normalizers b^P^^t or Y^x^oo^xbt to be negative. However, in our experiments, the latter was never negative and 
the former was very rarely negative; and, using real-valued observations with KDE (as in Section |3.5| l makes negative 
normalizers even less likely, since in this case the normalizer is a weighted sum of several estimated probabilities. In 
practice we recommend thresholding the normalizers with a small positive number, and not trusting probability estimates 
for a few steps if the normalizers fall below the threshold. 



5 



Note that the inference operations occur entirely in Mr. We mentioned earlier that parameterizing RR-HMM parameters as 
W x for all observations x casts it as a PSR of k dimensions. In fact the learning and inference algorithms for RR-HMMs 
proposed here have no dependence on the number of states m whatsoever, though other learning algorithms for RR-HMMs 
can depend on m (e.g. if they learn R and S directly). The RR-HMM formulation is intuitively appealing due to the idea 
of a large discrete state space with low-rank transitions, but this approach is also a provably consistent learning algorithm 
for PSRs in general, with finite-sample performance guarantees for the case where the PSR is an RR-HMM. Since PSRs 
are provably more expressive and compact than finite-state HMMs ]5]|T2], this indicates that we can learn a more powerful 
class of models than HMMs using this algorithm. 

3.3 Theoretical Guarantees 

The following finite sample bound on the estimated model generalizes analogous results from HKZ to the case of low -rank 
T. Theorem [2]bounds the L\ error in joint probability estimates from the learned model. This bound shows the consistency 
of the algorithm in learning a correct observable representation of the underlying RR-HMM, without ever needing to recover 
the high-dimensional parameters R, S, O of the latent representation. Note that our error bounds are independent of m, the 
number of hidden states; instead, they depend on k, the rank of the transition matrix, which can be much smaller than m. 
Since HKZ explicitly assumes a full-rank HMM transition matrix, and their bounds become vacuous otherwise, generalizing 
their framework involves relaxing this condition, generalizing the theoretical guarantees of HKZ and deriving proofs for 
these guarantees. 

Define tr^ (M) to denote the fc th largest singular value of matrix M, The sample complexity bounds depend polynomially 
on l/<7fc(P2,i) an d l/&k(OR). The larger er^p^i) is, the more well-separated are the dynamics from noise. The larger 
<7fc (OR) is, the more informative the observation is regarding state. For both these quantities, the larger the magnitude, the 
fewer samples we need to learn a good model. The bounds also depend on a term n Q (e), which is the minimum number 
of observations that account for (1 — e) of the total probability mass, i.e. the number of "important" observations. Recall 
that N is the number of independently sampled observation triples which comprise the training data, though as mentioned 
earlier we can also learn from a single long training sequence. 

The theorem holds under mild conditions. Some of these are the same as (or relaxations of) conditions in HKZ, namely 
that the prior ir is nonzero everywhere, and a number of matrices of interest (R, S, O, (U T OR)) are of rank at least k for 
invertibility reasons. The other conditions are unique to the low-rank setting, namely that S diag(7?)0 T has rank at least k, 
R has at least one column whose L2 norm is at most y/k/m, and the L\ norm of R is at most 1. The first of these conditions 
implies that the column space of S and the row space of O have some degree of overlap. The other two are satisfied, in the 
case of HMMs, by thinking of R as containing k linearly independent probability distributions along its columns (including 
a near-uniform column) and of S as containing the coefficients needed to obtain T from those columns. Alternatively, the 
conditions can be satisfied for an arbitrary R by scaling down entries of R and scaling up entries of S accordingly. However, 
this increases l/<Jk(OR), and hence we pay a price by increasing the number of samples needed to attain a particular error 
bound. See the Appendix (Section |7TT| ) for formal statements of these conditions. 

Theorem 2 [Generalization of HKZ Theorem 6] There exists a constant C > such that the following holds. Pick any 
< e, ?7 < 1 and t > 1. Assume the HMM obeys Conditions 3,4,5,6 and 7. Let e — <7fc(OP)<7fc(P2,i)e/ (4t\/k). Assume 

^ t 2 ( k k ■ n (e) \ , ... 

N ~ ° ' ? ' UiORWW + kWWWJ ' l0g(lA?) 
With probability > 1 — 77, the model returned by LEARNRR-HMMffc, N) satisfies 

^ I Pr[sGi, . . . ,x t ] -Pr[xi,...,^ t ]| < e 

Xl,.,.,Xt 

For the proof, see the Appendix (Section [74) >. 

3.4 Learning with Observation Sequences as Features 

The probability matrix Pi \ relates one past timestep to one future timestep, under the assumption that the vector of obser- 
vation probabilities at a single step is sufficient to disambiguate state (n > m and rank(O) = m). In system identification 
theory, this corresponds to assuming 1-step observability [9]. This assumption is unduly restrictive for many real-world 
dynamical systems of interest. More complex sufficient statistics of past and future may need to be modeled, such as the 
block Hankel matrix formulations for subspace methods |0[l7) to identify linear systems that are not 1-step observable. 
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For RR-HMMs, this corresponds to the case where n < m and/or rank(O) < m. Similar to the Hankel matrix formulation, 
we can stack multiple observation vectors such that each augmented observation comprises data from several, possibly 
consecutive, timesteps. The observations in the augmented observation vectors are assumed to be non-overlapping, i.e. all 
observations in the new observation vector at time t + 1 have larger time indices than observations in the new observation 
vector at time t. This corresponds to assuming past sequences and future sequences spanning multiple timesteps as events 
that characterize the dynamical system, causing Pi,P 2 .i and Pz,x,\ to be larger. Note that the x in P 3>x ,i still denotes a 
single observation, whereas the other indices in P\, P 2< i and P3, Xl i are now associated with events. For example, if we 
stack n consecutive observations, P 3fX equals the probability of seeing the i th n-length sequence, followed by the 
single observation x, followed by the j th n-length sequence. Empirically estimating this matrix consists of scanning for the 
appropriate subsequences i and j separated by observation symbol x, and normalizing to obtain the occurrence probability. 

P2,i and P3,x,i become larger matrices if we use a larger set of events in the past and future. However, stacking observations 
does not complicate the dynamics: it can be shown that the rank of P21 and P 3tX ^i cannot exceed k (see Section[8]in the 
Appendix for a proof sketch). Since our learning algorithm relies on an SVD of P 2 ,i, this means that augmenting the 
observations does not increase the rank of the HMM we are trying to recover. Also, since Pa^i is still an observation 
probability matrix with respect to a single unstacked observation x in the middle, the number of observable operators we 
need remains constant. Our complexity bounds successfully generalize to this case, since they only rely on Pi, P 2t i and 
P^,x,i being matrices of probabilities summing to 1 (for the former two) or to Pr[x2 = x] (for the latter), as they are here. 

The extension given above for learning HMMs with ambiguous observations differs from the approach suggested by HKZ, 
which simply substitutes observations with overlapping tuples of observations (e.g. P 2 ,i(j, i) — Pr[x 3 = j 2 ,x 2 = ji, %2 — 
*2> x i = There are two potential problems with the HKZ approach. First, the number of observable operators increases 
exponentially with the length of each tuple: there is one observable operator per tuple, instead of one per observation. 
Second, P 2 .i cannot be decomposed into a product of matrices that includes T, and consequently no longer has rank equal 
to the rank of the HMM being modeled. Thus, the learning algorithm could require much more data to recover a correct 
model if we use the HKZ approach. 



3.5 Learning with Real- Valued Observations 



The default RR-HMM formulation assumes discrete observations. However, since the model formulation converts the 
discrete observations into n-dimensional probability vectors, and the filtering, smoothing and learning algorithms we discuss 
all do the same, it is straightforward to model multivariate continuous data with Kernel Density Estimation [ 19|. 

This affects the learning algorithm and inference procedure as follows. Assume for ease of notation 
that the training data consists of N sets of three consecutive continuous observation vectors each, i.e., 
{(xi t i,^i t 2,xi t3 ), (x 2 , i,x 2 ,2,x 2 , 3), ■ ■ ■ , (xn,i,£n,2,x'n,3)}, though in practice we could be learning from a single long 
sequence (or several). Also assume for now that each observation vector contains a single raw observation, though this 
technique can easily be combined with the more sophisticated sequence-based learning and feature-based learning methods 
described above. Pick a kernel function K(-) and n kernel centers c.\ . . . c n . (In general we can use different kernels and 
centers for different feature vectors.) Let A be a bandwidth parameter that goes to zero at the appropriate rate in the limit. 

First compute n x 1 feature vectors (ifij)^-^, (ijjj)jLi, and (Cj)^Li> an d normalize each to sum to 1: 

[<f>j]i oc K(x h i - Ci) [ipj]i oc K(xj >2 - ?i) 

oc K(x Ji3 - Ci) [C 3 ]i oc K ((x Jt2 - Ci)/\) 

Note that for the second observation (in (Cj)yLi) we scale the kernel function by the bandwidth. Then, estimate the vector 
Pi and matrices P 2 i and P3 X: i (for x = &i, . . . ,c n ) from data: 

N N 

j=i 3=1 



1 N _ 

For x = ci, . . . , c„: P 3tXA = — ^[QxZj^j 



We compute n 'base' observable operators B Cl , . . . , B Cri from the estimated probability matrices, as well as vectors bi and 



boo, using algorithm Learn-RR-HMM (Section 3.1 1. Given these parameters, filtering for a sequence (xi, . . . , x T ) now 
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Figure 2: Learning discrete RR-HMMs. The three figures depict the actual eigenvalues of three different RR-HMM tran- 
sition matrices, and the eigenvalues (95% error bars) of the sum of RR-HMM observable operators estimated with 10, 000 
and 100, 000 training observations. (A) A 3-state, 3-observation, rank 2 RR-HMM. (B) A full-rank, 3-state, 2-observation 
HMM. (C) A 4-state, 2-observation, rank 3 RR-HMM. 



proceeds as follows: 

For t = 1, . . . , r: 

Compute and normalize [a t ]i oc K ((xt — Cj)/A) . 

n 

h B ^ 

oi+i — 



b 00 B at b t 

Our theoretical results carry over to the KDE case with modifications described in the RR-HMM document. Essentially, the 
bound still holds for predicting functions of <7i, <T2, . . . , &u though we do not yet have results connecting this bound to the 
error in estimating probabilities of raw observations. 



4 Experimental Results 

We designed several experiments to evaluate the properties of RR-HMMs and the learning algorithm both on synthetic 
and on real-world data. The first set of experiments (Section |4. \\ tests the ability of the spectral learning algorithm to 



recover the correct RR-HMM. The second experiment (Section 4.2 1 evaluates the representational capacity of the RR- 
HMM by learning a model of a video that requires both competitive inhibition and smooth state evolution. The third set of 
experiments (Section |4~3] > tests the model's ability to leam, filter, predict, and simulate video captured from a robot moving 
in an indoor office environment. 



4.1 Learning Synthetic RR-HMMs 

First we evaluate the unbiasedness of the spectral learning algorithm for RR-HMMs on 3 synthetic examples. In each case, 
we build an RR-HMM, sample observations from the model, and estimate the model with the spectral learning algorithm 
described in Section[3] We compare the eigenvalues of B — ^2 X B x in the learned model to the eigenvalues of the transition 
matrix T of the true model. B is a similarity transform of S ■ R which therefore has the same non-zero eigenvalues as 
T = RS, so we expect the estimated eigenvalues to converge to the true eigenvalues with enough data. This is a necessary 
condition for unbiasedness but not a sufficient one. See Section|9]in Appendix for parameters of HMMs used in the examples 
below. 

Example 1: An RR-HMM We examine an HMM with m = 3 hidden states, n = 3 observations, a full-rank observation 
matrix and a k = 2 rank transition matrix. Figure [2jA) plots the true and estimated eigenvalues for increasing size of 
dataset, along with error bars, suggesting that we recover the true dynamic model. 
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Figure 3: The clock video texture simulated by a HMM, a stable LDS, and a RR-HMM. (A) The clock modeled by a 
10-state HMM. The manifold consists of the top 3 principal components of predicted observations during simulation. The 
generated frames are coherent but motion in the video is jerky. (B) The clock modeled by a 10-dimensional LDS. The 
manifold indicates the trajectory of the model in state space during simulation. Motion in the video is smooth but frames 
degenerate to superpositions. (C) The clock modeled by a rank 10 RR-HMM. The manifold consists of the trajectory of the 
model in the low dimensional subspace of the state space during simulation. Both the motion and the frames are correct. 



Example 2: A 2-step-Observable HMM We examine an HMM with m — 3 hidden states, n = 2 observations, and a 
full-rank transition matrix (see Appendix for parameters). This HMM violates the m < n condition. The parameters of 
this HMM cannot be estimated with the original learning algorithm, since a single observation does not provide enough 



information to disambiguate state. By stacking 2 consecutive observations (see Section 3.4 1, however, the spectral learning 
algorithm can be applied successfully (Figure|2|B)). 

Example 3: A 2-step-Observable RR-HMM We examine an HMM with m = 4 hidden states, n = 2 observations, 
and a k = 3 rank transition matrix (see Appendix for parameters). In this example, the HMM is low rank and multiple 
observations are required to disambiguate state. Again, stacking two consecutive observations in conjunction with the 
spectral learning algorithm is enough to recover good RR-HMM parameter estimates (Figure[2|C)). 



4.2 Competitive Inhibition and Smooth State Evolution in Video 

We model a clock pendulum video consisting of 55 frames (with a period of ~ 22 frames) as a 10-state HMM, a 10- 
dimensional LDS, and a rank 10 RR-HMM with 4 stacked observations. Note that we could easily learn models with 
more than 10 latent states/dimensions; we limited the dimensionality in order to demonstrate the relative expressive power 
of the different models. For the HMM, we convert the continuous data to discrete observations by 1-NN on 25 kernel 
centers sampled sequentially from the training data. We trained the resulting discrete HMM using EM. We learned the LDS 
directly from the video using subspace ID with stability constraints |20| using a Hankel matrix of 10 stacked observations. 
We trained the RR-HMM by stacking 4 observations, choosing an approximate rank of 10 dimensions, and learning 25 
observable operators corresponding to 25 Gaussian kernel centers. We simulate a series of 500 observations from the model 
and compare the manifolds underlying the simulated observations and frames from the simulated videos (Figure [3). The 
small number of states in the HMM is not sufficient to capture the smooth evolution of the clock: the simulated video is 
characterized by realistic looking frames, but exhibits jerky irregular motion. For the LDS, although the 10-dimensional 
subspace captures smooth evolution of the simulated video, the system quickly degenerates and individual frames of video 
are modeled poorly (resulting in superpositions of pendulums in generated frames). For the RR-HMM, the simulated video 
benefits from both smooth state evolution and competitive inhibition. The manifold in the 10-dimensional subspace is 
smooth and structured and the video is realistic. The results demonstrate that the RR-HMM has the benefits of smooth state 
evolution and compact state space of a LDS and the benefit of competitive inhibition of a HMM. 



4.3 Filtering, Prediction, and Simulation with Robot Vision Data 

We compare HMMs, LDSs, and RR-HMMs on the problem of modeling video data from a mobile robot in an indoor 
environment. A video of 2000 frames was collected at 6 Hz from a Point Grey Bumblebee2 stereo camera mounted on a 
Botrics Obot dlOO mobile robot platform circling a stationary obstacle (Figure[4]A)) and 1500 frames were used as training 



9 



A. 



Example Images 




















B. xl ° 

8 



•2 6 

O 

^ 5 

CD 

S-H 

$ 3 



L 

1 






■ — ■ Mean 

— — Last 

RR-HMM 

LDS 

"»HMM 



Environment 



10 20 30 40 50 60 70 80 90 100 

Prediction Horizon 



Figure 4: (A) Sample images from the robot's camera. The figure below depicts the hallway environment with a central 
obstacle (black) and the path that the robot took through the environment while collecting data (the red counter-clockwise 
ellipse) (B) Squared error for prediction (1, . . . , 100 steps out in future) with different estimated models and baselines, 
averaged over different initial filtering durations (1, . . . , 250). 



data for each model. Each frame from the training data was reduced to 100 dimensions via SVD on single observations. 
Using this training data, we trained an RR-HMM (k = 50, n = 1500) using spectral learning with sequences of 20 
continuous observations (Section [3. 4) and KDE with Gaussian kernels (Section [33} with 1500 centers, a 50-dimensional 
LDS using Subspace ID with Hankel matrices of 20 timesteps, and a 50-state HMM with 1500 discrete observations using 
EM run until convergence. For each model, we performed filtering for different extents t± = 100, 101, . . . ,250, then 
predicted an image which was a further t-x steps in the future, for = 1,2 ... , 100. The squared error of this prediction 
in pixel space was recorded, and averaged over all the different filtering extents t\ to obtain means which are plotted in 
Figure |4jB). As baselines, we plot the error obtained by using the mean of filtered data as a predictor ('Mean'), and the 
error obtained by using the last filtered observation ('Last'). 

Both baselines perform worse than any of the more complex algorithms (though as expected, the 'Last' predictor is a good 
one-step predictor), indicating that this is a nontrivial prediction problem. The LDS does well initially (due to smoothness), 
and the HMM does well in the longer run (due to competitive inhibition), while the RR-HMM performs as well or better at 
both time scales since it models both the smooth state evolution and competitive inhibition in its predictive distribution. In 
particular, the RR-HMM yields significantly lower prediction error consistently for the duration of the prediction horizon 
(100 timesteps, i.e. 16 seconds). 



5 Related Work 

5.1 Predictive State Representations 

Predictive State Representations (PSRs) llUED an d Observable Operator Models (OOMs) |5| model sequence probabil- 
ities as a product of observable operator matrices. This idea, as well as the idea of learning such models using linear 
algebra techniques, originates in the literature on multiplicity automata and weighted automata QIHllSl- Despite recent 
improvements 12T1I221 . practical learning algorithms for PSRs and OOMs have been lacking. RR-HMMs and its spectral 
learning algorithm are also closely related to methods in subspace identification ll9l [T7ll in control systems for learning LDS 
parameters, which use SVD to determine the relationship between hidden states and observations. 

As pointed out earlier, the spectral learning algorithm presented here learns PSRs. We briefly discuss other algorithms 
for learning PSRs from data. Several learning algorithms for PSRs have been proposed ll23l [T21 l24l . It is easier for PSR 
learning algorithms to return consistent parameter estimates because the parameters are based on observable quantities. 
ifTTI develops an SVD-based method for finding a low-dimensional variant of PSRs, called Transformed PSRs (TPSRs). 
Instead of tracking the probabilities of a small number of tests, TPSRs track a small number of linear combinations of a 
larger number of tests. This allows more compact representations, as well as dimensionality selection based on examining 
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the singular values of the decomposed matrix, as in subspace identification methods. Note that nonlinearity can be encoded 
into the design of core tests. [25 1 introduced the concept of e-tests in PSRs that are indicator functions of aggregate 
sets of future outcomes, e.g. all sequence of observations in the immediate future that end with a particular observation 
after k timesteps. In general, tests in discrete PSRs can be indicator functions of arbitrary statistics of future events, thus 
encoding nonlinearities that might be essential for modeling some dynamical systems. Recently, Exponential Family PSRs 
(EFPSRs) [21 1 were introduced as an attempt to generalize the PLG model to allow general exponential family distributions 
over the next N observations. In the EFPSR, state is represented by modeling the parameters of a time-varying exponential 
family distribution over the next N timesteps. This allows graphical structure to be encoded in the distribution, by choosing 
the parameters accordingly. The justification for choosing an exponential family comes from maximum entropy modeling. 
Though inference and parameter learning are difficult in graphical models of non-trivial structure, approximate inference 
methods can be utilized to make these problems tractable. Like PLGs, the dynamical component of EFPSRs is modeled by 
extending and conditioning the distribution over time. However, the method presented [21] has some drawbacks, e.g. the 
extend-and-condition method is inconsistent with respect to marginals over individual timesteps between the extended and 
un-extended distributions. 

5.2 Hybrid Models, Mixture Models and other recent approaches 

RR-HMMs and their algorithms are also related to other hybrid models. Note that previous models of the same name 
(e.g. [26 1) address a completely different problem, i.e. reducing the rank of the Gaussian observation parameters. Since 
shortly after the advent of LDSs, there have been attempts to combine the discrete states of HMMs with the smooth dy- 
namics of LDSs. We perform a brief review of the literature on hybrid models; see G71 for a more thorough review. Il28ll 
formulates a switching LDS variant where both the state and observation variable noise models are mixture of Gaussians 
with the mixture switching variable evolving according to Markovian dynamics, and derives the (intractable) optimal filter- 
ing equations where the number of Gaussians needed to represent the belief increases exponentially over time. They also 
propose an approximate filtering algorithm for this model based on a single Gaussian. |29| proposes learning algorithms 
for an LDS with switching observation matrices. 1 30 1 reviews models where both the observations and state variable switch 
according to a discrete variable with Markov transitions. Hidden Filter HMMs (HFHMMs) ||3D combine discrete and real- 
valued state variables and outputs that depend on both. The real-valued state is deterministically dependent on previous 
observations in a known manner, and only the discrete variable is hidden. This allows exact inference in this model to be 
tractable. [32| formulates the Mixture Kalman Filter (MKF) model along with a filtering algorithm, similar to [28| except 
that the filtering algorithm is based on sequential Monte-Carlo sampling. 

The commonly used HMMs with mixture-model observations (e.g., Gaussian mixture) are a special case of RR-HMMs. A k- 
state HMM where each state corresponds to a Gaussian mixture of m observation models of n dimensions each is subsumed 
by a fc-rank RR-HMM with m distinct continuous observations of n dimensions each, since the former is constrained to 
be non-negative and < 1 in various places (the fc-dimensional transition matrix, the fc-dimensional belief vector, the matrix 
which transforms this belief to observation probabilities) where the latter is not. 

Switching State-Space Models (SSSMs) [27 1 posit the existence of several real-valued hidden state variables that evolve 
linearly, with a single Markovian discrete-valued switching variable selecting the state which explains the real-valued ob- 
servation at every timestep. Since exact inference and learning are intractable in this model, the authors derive a structured 
variational approximation that decouples the state space and switching variable chains, effectively resulting in Kalman 
smoothing on the state space variables and HMM forward-backward on the switching variable. In their experiments, the 
authors find SSSMs to perform better than regular LDSs on a physiological data modeling task with multiple distinct un- 
derlying dynamical models. HMMs performed comparably well in terms of log-likelihood, indicating their ability to model 
nonlinear dynamics though the resulting model was less interpretable than the best SSSM. More recently, models for non- 
linear time series modeling such as Gaussian Process Dynamical Models have been proposed [33 1] - However, the parameter 
learning algorithm is only locally optimal, and exact inference and simulation are very expensive, requiring MCMC over a 
long sequence of frames all at once. This necessitates the use of heuristics for both inference and learning. Another recent 
nonlinear dynamic model is ll34l . which differs greatly from other methods in that it treats each component of the dynamic 
model learning problem separately using supervised learning algorithms, and proves consistency on the aggregate result 
under certain strong assumptions. 

6 Discussion 

The spectral learning algorithm blurs the line between latent variable models and PSRs. PSRs were developed with a focus 
on the problem of an agent planning actions in a partially observable environment. More generally, there are many scenarios 
in sequential data modeling where the underlying dynamical system has inputs. The inference task for a learned model is 
then to track the belief state while conditioning on observations and incorporating the inputs. The input-output HMM (10- 
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HMM) 11351 is a conditional probabilistic model which has these properties. A natural generalization of this work is to the 
task of learning RR-HMMs with inputs, or controlled PSRs. We recently carried out this generalization to controlled PSRs; 
details can be found in iflOl . 

The question of proving containment or equivalence of RR-HMMs with respect to PSRs is of theoretical interest. The 
observable representation of an RR-HMM is a Transformed PSR (TPSR) [11], so every RR-HMM is a PSR; it remains to 
be seen whether every PSR corresponds to some RR-HMM (possibly with an infinite number of discrete hidden states) as 
well. The idea that "difficult" PSRs should somehow correspond to RR-HMMs with very large or infinite state space is 
intuitively appealing but not straightforward to prove. Another interesting direction would be to bound the performance of 
the learning algorithm when the underlying model is only approximately a reduced-rank HMM, much as the HKZ algorithm 
includes bounds when the underlying model is approximately an HMM (l |. This would be useful since in practice it is more 
realistic to expect any underlying system to not comply with the exact model assumptions. 

The positive realization problem, i.e. obtaining stochastic transition and observation matrices from the RR-HMM observable 
representation, is also significant, though the observable representation allows us to carry out all possible HMM operations. 
HKZ describes a method based on [36] which, however, is highly erratic in practice. In the RR-HMM case, we have the 
additional challenge of firstly computing the minimal rn for which a positive realization exists, and since the algorithm 
learns PSRs there is no guarantee that a particular set of learned parameters conforms exactly to any RR-HMM. On the 
applications side, it would be interesting to compare RR-HMMs with other dynamical models on classification tasks, as 
well as on learning models of difficult video modeling and graphics problems for simulation purposes. More elaborate 
choices of features may be useful in such applications, as would be the usage of high-dimensional or infinite-dimensional 
features via Reducing Kernel Hilbert Spaces (RKHS). 
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7 Appendix I: Proofs 



The proof of Theorem [2]relies on Lemmas 18 and 24 We start off with some preliminary results and build up to proving 
the main theorem and its lemmas below. 

A remark on norms : The notation \\X\\ p for matrices X 6 jj mx » l denotes the operator norm max nSr^ for vector u^O. 
Specifically, ||X|| 2 denotes L2 matrix norm (also known as spectral norm), which corresponds to the largest singular value 

( m n X 1 / 2 

o\{X). Frobenius norm is denoted by = [YmLi Y^j=i -^Hj ) ■ The notation WX^ for matrices denotes the L\ 

matrix norm which corresponds to maximum absolute column sum max c 2~2iLi \Xi C \- The definition of ||a;|| for vectors 
x e K ra is the standard distance measure (X)"=i ^f) 1 ^- 

7.1 Preliminaries 

The following conditions are assumed by the main theorems and algorithms. 

Condition 3 [Modification ofHKZ Condition 1] ir > element-wise, T has rank k (i.e. R and S both have rank k) and O 
has rank at least k. 

The following two conditions on R can always be satisfied by scaling down entries in R and scaling up S accordingly. 
However we want entries in R to be as large as possible under the two conditions below, so that ak(U T OR) is large and 
l/ak(U T OR) is small to make the error bound as tight as possible (Theorem Hence we pay for scaling down R by 
loosening the error bound we obtain for a given number of training samples. 

Condition 4 | \R\ | x < 1. 

Condition 5 For some column 1 < c < k of R, it is the case that \\R\-, c]|L < \Jkjra. 

The above two conditions on R ensure the bounds go through largely unchanged from HKZ aside from the improvement 
due to low rank k. The first condition can be satisfied in a variety of ways without loss of generality, e.g. by choosing 
the columns of R to be any k independent columns of T, and S to be the coefficients needed to reconstruct T from R. 
Intuitively, the first condition implies that R does not overly magnify the magnitude of vectors it multiplies with. The 
second one implies a certain degree of uniformity in at least one of the columns of R. For example, the uniform distribution 
in a column of R would satisfy the constraint, whereas a column of the identity matrix would not. This does not imply that 
T must have a similarly near-uniform column. We can form R from the uniform distribution along with some independent 
columns of T. 

The observable representation depends on a matrix U £ K nxfe that obeys the following condition: 
Condition 6 [Modification of HKZ Condition 2] U T OR is invertible. 

This is analogous to the HKZ invertibility condition on U t O, since OR is the matrix that yields observation probabilities 
from a low-dimensional state vector. Hence, U defines a fc-dimensional subspace that preserves the low-dimensional state 
dynamics regardless of the number of states m. 

Condition 7 Assume that S diag(7?)0 T has full row rank (i.e. k). 

This condition amounts to ensuring that that the ranges S and O, which are both at least rank k, overlap enough to preserve 
the dynamics when mapping down to the low-dimensional state. As in HKZ, the left singular vectors of P 2 .i gi ye us a valid 
U matrix. 

Lemma 8 [Modification of HKZ Lemma 2] Assume Conditions^and^ Then, rank(P2.i) = k. Also, ifU is the matrix of 
left singular vectors of P2,x corresponding to non-zero singular values, then range (U) = range(Oi?), so U € R nxk obeys 
Condition® 
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Proof: From its definition, we can show P^i can be written as a low-rank product of RR-HMM parameters: 

[P2,i]i,j = Pr[x 2 = i,xi = j] 

m m 

= Pr[x*2 = i, $x = j, hi = a,h\ = 6] (marginalizing hidden states /i) 

a=l 6=1 
m m 

= ^^Pr[x 2 = i|/i 2 = a] Pr[/i 2 = a|/ii = 6]Pr[.Ti = j\hi = b] Pr[/ii = 6] 

0=1 6=1 

m m 

= ^ ^ O ia T ab 7f b [0 T ] &:) 

a=l 6=1 

=*> P 2 .i = OT diag(7?)0 T 

= ORS diag(7r)0 T (7.1) 
Thus range(P2,i) Q range(OP). This shows that rank(P2,i) < rank(OP). 

By Condition^ 5 diag(7r)0 T has full row rank, thus S diag(7r)0 T (5 diag(7?)0 T ) + = Ikxk- Therefore, 

OR = P 2 ,i{S diag(7?)0 T )+ (7.2) 
which implies range(OP) C range(P2,i), which in turn implies rank(OP) < rank(P2,i). 

Together this proves that rank(P2 j i) = rank(OP), which we can show to be k as follows: Condition [3] implies that 
rank (U 1 OR) = k, and hence rank(OP) > k. Since OR E R mxk , rank(OP) < k. Therefore, rank(OP) = k. Hence 
rank(P2 ! i) = k. 

Since range([7) = range(P2.i) by definition of singular vectors, this implies range(P) — range(OP). Therefore, 
(U T OR) is invertible and hence U obeys Condition]^ □ 

The following lemma shows that the observable representation {boo, bi, Pi, . . . , B n } is linearly related to the true HMM 
parameters, and can compute the probability of any sequence of observations. 

Lemma 9 I Modification of HKZ Lemma 3] (Observable HMM Representation). Assume Condition^on the RR-HMM and 
Conditionpjon the matrix U G R" xfe . Then, the observable representation of the RR-HMM (Definition 1 of the paper) has 
the following properties: 

1. bi = (U t OR)tti = (U t O)tt, 

2. bZ = lJ n R(U T ORr\ 

3. Forallx=l,...,n: B x = (U T OR)W x (U T OR)- 1 

4. For any timet: Vx\x\j\ — fy^P^i&i 
Proof: 

1 . We can write p as Oir, since 

= Pr[zi = i] 

m 

= Pr[a;i = i\h\ = a] Pr[/ii = a] 

a=l 
m 

a=l 

Combined with the fact that bi = U T Pi by definition, this proves the first claim. 

2. Firstly note that Pj = T^T diag(7?)0 T , since 

PJ = ^o J 

= P m diag(7?)0 T 

= ZT diag(7r)0 T (since P m T = Tjj 
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This allows us to write Pi in the following form: 

PJ = F m T diag(^)0 T 

= T^RS diag(7r)0 T 

= l T m R(U T OR)- 1 (U T OR)S diag(7r)0 T 

= T^ n R{U T OR)- 1 U T P 2 ,i (by equation $71) ) 

By equation ((TT}, U 1 ' P 2 .i = (U T OR)S diag(7?)0 T . Since (U T OR) is invertible by Condition |6j and 
S diag(7?)0 T has full row rank by Condition]?] we know that (U T P 2t i) + exists and 

U T P 2A (U T P 2A ) + = I kxk (7.3) 

Therefore, 

bl = P 1 T (C/ T P 2 , 1 )+ = lj n R(U T OR)- 1 (U T P 2 ,i)(U T P 2 ,i) + = llR^OR)- 1 
hence proving the second claim. 
3. The third claim can be proven by first expressing P 3 , x .i as a product of RR-HMM parameters: 

[P3,x,l]ij = Pr[x 3 =i,X2 = x, xi = j] 



m m m 



^ X] X Pl '\- X 3 = h x 2= x, Xi = j, h 3 = a,h 2 = b, hi = 



a=l 6=1 c=l 
m m tci 



a=l 6=1 c=l 

Pr[ft, 2 = b\hi = c] Pr[/ii = c] Pr[:ri = = c] 



m m m 



= EEE O m [^J ab T bc 7f c [0 T ] CJ 

a=l 6=1 c=l 

P 3)X ,i = OA x T diag(^)0 T 
This can be transformed as follows: 

P 3 ,x,i = OA x RS diag(7?)0 T 

= OA x R(U T OR)- 1 (U T OR)S diag(7r)0 T 

= OA x R(U T OR)- 1 U T (OT diag(7r)0 T ) 

= OA x R(U T OR)- 1 U T P 2 , 1 (by equation (71) ) 

and then plugging in this expression into the definition of B x , we obtain the required result: 

B x = (U T P 3 , xA )(U T P 2A )+ (by definition) 
= (U T 0)A x R(U T OR)-\U T P 2 , 1 )(U T P 2 , 1 ) + 
= (U T 0)A x R(U T OR)- 1 (by equation {73])) 
= (U T OR) (S di&g(O x ,)R) (l/TOR)- 1 
= (U T OR)W x (U T OR)- 1 



4. Using the above three results, the fourth claim follows from equation 1 in Section 2 in the paper: 
Pr[a;i, ...,x t ] 

= r v m RW xt ...w xl n 

= lJ n R(U T OR)- 1 (U T OR)W Xt (U J OR)- 1 (U T OR)W Xt _ 1 (U T OR)- 1 . . . 

. . . (U t OR)W Xi (U t OR)- 1 (U t OR)tti 
= bL B xt ■ ■ ■ B Xl bi 
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□ 

In addition to b\ above, we define normalized conditional 'internal states' b t that help us compute conditional probabilities. 
These internal states are not probabilities. In contrast to HKZ where these internal states are m-dimensional vectors, in our 
case the internal states are fc-dimensional i.e. they correspond to the rank of the HMM. As shown above in Lemma|9] 

h = (U T OR)TTt = (C/ T 0)vr 
In addition for any t > 1, given observations Xi-.t-i with non-zero probability, the internal state is defined as: 

b t = bt{x lxt - x ) = (7.4) 

For t = 1 the formula is still consistent since b^bi = Tj n R(U T ORy 1 (U T OR)TTi = l^R^i = T^n = 1. 
Recall that HMM and RR-HMM parameters can be used to calculate joint probabilities as follows: 

Pr[xi, x t ] = l y m A Xt A Xt _ 1 ■ ■ ■ A Xl n 

= l^RS dia,g(0 Xt ,)RS dia,g(0 Xt _ u .)R---S diag(O ai> .)7r 

= F m R(S diag(0 Xti .) J R) (S diag(0 Xt _ u .)R)---(S diag(0 Xl ,.)R) Sir 

= T^RW Xt ...W X1 n (by definition of W x , n) (7.5) 

The following Lemma shows that the conditional internal states are linearly related, and also shows how we can use them 
to compute conditional probabilities. 

Lemma 10 [Modification of HKZ Lemma 4] (Conditional Internal States) Assume the conditions of Lemma^hold, i.e. 
Conditions^and^hold. Then, for any time t: 

1. (Recursive update) If Pr[xi, . . . , Xt] > 0, then 

r B xt b t 

°t+l — — 

blB Xt b t 

2. (Relation to hidden states) 

b t = (U r OR)k(x v . t -i) = (U 1 '0)ht(x 1:t -i) 

where [/i^ (aci:t— 1)]< = Pr[/it = «|^i:t-i] ' s defined as the conditional probability of the hidden state at time t given 
observations x\-t~\, and l t (xi-,t~i) is its low-dimensional projection such that h t {x\;t~i) — Rlt(xi-,t-i). 

3. ( Conditional observation probabilities) 

Pi[x t \xi :t -i] = b^B^bt 
Proof: The first proof is direct, the second follows by induction. 



1 . The t = 2 case b 2 = Baj i fcl ^ [ s true by definition (equation 



7.4 1. For t > 3, again by definition of fe t+1 we have 



b 



t+i 



B Xt B Xt _ 1:1 bi 

i ui 

B b 

— — (by equation %Ffy) 

B b 

_ Xt ^ (by equation (74)) 
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2,3. The base case for claim 2 holds by Lemma[9] since h\ = 7?, h = -Rt? and 6i = (U T OR)if. For claim 3, the base 
case holds since b^Bx^i = Tl n RW Xl iri by Lemma|9] which equals Pr[a;i] by equation ( |7.5| >. The inductive step 



is: 



bt+i — - — (by claim 1 above) 

oLB Xt b t 

B Xt (U T OR)k t , . , 

— (by inductive hypothesis) 



Pr[x t \xi:t-i} 
= (U T OR)W Xt l t 

Pr[x t \xi:t-i] 
= (U T Q)A Xt h t 
Pr[x t \xi;t-i] 

Now by definition of A Xt h t , 



(by Lemma [9]) 
(••■ RW Xt l t = RS ding(0 Xt7 .)Rl t = Ajt) 



Pr[xt\xi :t -i] 

Pi[x t \x 1: t-i] 
= (U T 0)h t+1 (x 1:t ) 
= (U T OR% +1 (x 1:t ) 
This proves claim 2, using which we can complete the proof for claim 3: 

blB Xt+1 b t+1 = ]J n R(U T OR)- 1 (U T OR)W Xt (U T OR)-% +1 (by Lemma[9) 
= tj n RW Xt (U T OR)- 1 (U T OR)l t+ i (by claim 2 above) 

= r m RS &ia,g{0 Xt ,)Rl t+l (by definition of W Xt ) 
= l^T diag(<9 Xt ,.)£ t+1 
— ^ln.A Xt h t+ i 
Again by definition of A Xt h t +\, 

m rn 

^oo^t+i^t+i = }^} j Pi[x t+1 \h t+ i = a]Pr[/i t+ i = aj/it = 6] Pr[/i t = 6|xi :t ] 
a=i 6=1 

m m 

= Pr[2^ + i, /i t+ i = a, = 6|aci: t ] 

a=l 6=1 



□ 



Remark 11 If U is the matrix of left singular vectors of Pi\ corresponding to non-zero singular values, then U is the 
observable-representation analogue of the observation probability matrix O in the sense that, given a conditional state b t , 
Pr[x t = i\%i : t—i] = [Ub t ]i in the same way as Pr[x t = = [Oh t ]ifor a conditional hidden state h t . 

Proof: Since range (L/) = range(Oi?) (Lemma 2), and UU T is a projection operator to range({7), we have UU T OR = 
OR, so Ub t = U(U T OR)l t = ORl t = Oh t . □ 

7.2 Matrix Perturbation Theory 

We take a diversion to matrix perturbation theory and state some standard theorems from Steward and Sun (1990) ll37ll 
and Wedin (1972) ll38l which we will use, and also prove a result from these theorems. The following lemma bounds the 
L2-norm difference between the pseudoinverse of a matrix and the pseudoinverse of its perturbation. 
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Lemma 12 (Theorem 3.8 of Stewart and Sun (1990) 071/ ) Let A G 

1 



1 ', with m > n, and let A = A + E. Then, 



A+-A^ 



v 5 r 1 1 .1 1 r 2 

< • max < L4 T 

2 1,2 



A+ 



\E\\ 



The following lemma bounds the absolute differences between the singular values of a matrix and its perturbation. 

Lemma 13 (Theorem 4.11 of Stewart and Sun (1990) 071/ ). Let A G R mxn with m > n, and let A = A 
singular values of A and A are (crj > . . . > a n ) and (a% > . . . > a n ), respectively, then 



E. If the 



< \\E\L i = l, 



,n 



Before the next lemma we must define the notion of canonical angles between two subspaces: 

Definition 14 (Adapted from definition 4.35 of Stewart (1998) [39]) Let X and Y be matrices whose columns comprise 
orthonormal bases of two p-dimensional subspaces X and y respectively. Let the singular values of X T Y ( where X T 
denotes the conjugate transpose, or Hermitian, of matrix X) be 71, 72, ... , 7 p . Then the canonical angles 9i between X and 
y are defined by 

6i = cos -1 7,:, i=l,2,...,p 
The matrix of canonical angles O is defined as 

e(x,y) = diag(0i,0a,...,0 p ) 

Note that Vz 7$ G [0, 1] in the above definition, since 71 (assuming it's the highest singular value) is no greater than 
a\(X T )ai(Y) < 1 • 1 = 1, and hence cos -1 7$ is always well-defined. 

For any matrix A, define A± to be the orthogonal complement of the subspace spanned by the columns of A. For example, 
any subset of left singular vectors of a matrix comprise the orthogonal complement of the matrix composed of the remaining 
left singular vectors. The following lemma gives us a convenient way of calculating the sines of the canonical angles between 
two subspaces using orthogonal complements: 

Lemma 15 (Theorem 4.37 of Stewart (1998) 091/ ) Let X and Y be n x p matrices, with n > p, whose columns comprise 
orthonormal bases of two p-dimensional subspaces X and y respectively. Assume Xj_, Y± G R" x "" p such that [X Xj_] 
and [Y Yj_] are orthogonal matrices. The singular values ofY^X are the sines of the canonical angles between X and y. 

The following lemma bounds the L 2 -norm difference between the sine of the canonical angle matrices of the range of a 
matrix and its perturbation. 



Lemma 16 ( KM . Theorem 4.4 of Stewart and Sun (1990) [37]). Let A e 
decomposition (Ui, U%, U3, £1, £2, Vj., V2): 



with m > n, with the singular value 



UJ 



A [ Vi V 2 ] = 



£1 







£2 





Let A = A + E, with analogous SVD (U\, U2, U3, Si, £2, S3, Vj., Vi)- Let be the matrix of canonical angles between 
range([/i) and range(Z7i), and be the matrix of canonical angles between range(Vi) and range^i). If there exists 
8 > 0, a > such that min cr(Si) > a + 5 and max<7(S 2 ) < a > then 

\\E\\ 

max{||sin$|| 2 , ||sin6|| 2 } < 



The above two lemmas can be adapted to prove that Corollary 22 of HKZ holds for the low-rank case as well, assuming 
that the perturbation is bounded by a number less than Uk- The following lemma shows that (1) the k th singular value of 
a matrix and its perturbation are close to each other, and (2) that the subspace spanned by the first k singular vectors of a 



matrix is nearly orthogonal to the subspace spanned by the (k + 1) 
matrix product of their bases being bounded. 



ih 



, m singular vectors of its perturbation, with the 



Corollary 17 [Modification of HKZ Corollary 22] Let A G R mxn , with m > n, have rank k < n, and let U G W nxk be 
the matrix ofk left singular vectors corresponding to the non-zero singular values o\ > . . . > Ok > of A. Let A = A+E. 
Let U G M. mxk be the matrix ofk left singular vectors corresponding to the largest k singular values U\ > . . . > Uf. of A, 
and let U± G Jj mx ( m - fc ) }j e the remaining left singular vectors. Assume \\E\\ 2 < eak for some e < 1. Then: 
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1- <?*>(!- e)a k . 



UjU 



< \\E\\ 2 /cr k . 



2. 

Proof: 

1. FromLemmafT3 



\a k -a k \ < \\E\\ 2 
\cfk — cr fc | < ecr fe 
Ofc - Cfc > -ecr fc 

5fe > (1 - e)a k 



which proves the first claim. 



2. Recall that by Lemma 15 if $ is a matrix of all canonical angles between range(P2,i) and range(P2,i), then sin $ 
contains all the singular values of UjU along its diagonal. 
Also recall that the L 2 norm of a matrix is its top singular value. Then, 

||sin$|| 2 = <7i(sin<l>) (by definition) 

= max diag(sin $) (since sin $ is a diagonal matrix) 

= o-i {UjU) (by LemmafB i 



UjU 



(by definition) 



Invoking Lemma 
with ||sin$|U = 



16 



with the parameter values 8 = a k and a = yields |jsin<I>|| 2 < \\E\\ 2 ja k . Combining this 



(UlU) 



proves claim 2. 



□ 



7.3 Supporting Lemmas 

In this section we develop the main supporting lemmas that help us prove Theorem [2] 
7.3.1 Estimation Errors 

We define £i,£2,i and e-$ :X ,i as sampling errors for P\,P-z t i and Pa^i respectively: 



ei = 

£2,1 = 
£3.x, 1 == 



Pi -Pi 
?2,1 - P 



2.1 



P- 



3,a;,l 



P 



3,z,l 



for x 



(7.6a) 
(7.6b) 
(7.6c) 



Lemma 18 [Modification of HKZ Lemma 8] If the algorithm independently samples N observation triples from the HMM, 
then with probability at least 1 — r\: 



ei < 
£2,1 < 
maxe 3iIi i < 





In 


3 






N 




n 






3 




In 




N 




n 






3 




In 




N 




n 



y^£3,a,l < min ^ 



N 
N 
N 
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Before proving this lemma, we need some definitions and a preliminary result. First, we restate McDiarmid's Inequal- 
ity IH: 



Theorem 19 Let Z\ , . . . , Z rn be independent random variables all taking values in the set Z. Let Ci be some positive real 
numbers. Further, let f : Z m i— > R be a function of Z\, . . . , Z m f/jaf satisfies Vi, Vzi, . . . , z m , z[ £ Z, 



|/(zi,...,z l; ...,z m )-/(zi, 



, Zj , . . . , z m ) I < Cj 



Then for all e > 0, 



l'v[/ - '-£[/] .■<] :exp( =J^ 2 . 



Assume z is a discrete random variable that takes on values in 1, . . . , d. The goal is to estimate the vector q = [Pr(z = 
from TV i.i.d. samples Zi (i = 1, . . . , N). Let e 3 denote the j th column of the d x d identity matrix. For i — 1, . . . , N, 
suppose q. L is a column of the d x d identity matrix such that — e Zi . In other words, the zf l component of qi is 1 and 
the rest are 0. Then the empirical estimate of (fin terms of qi is q — J2i=i Qi/N. 

Each part of Lemma [T~8] corresponds to bounding, for some q, the quantity 

\\q-q\\l ■ 



We first state a result based on McDiarmid's inequality (Theorem[T9j): 

Proposition 20 [Modification of HKZ Proposition 19] For all e > and q, q and N as defined above: 



Pr ( ||g — g|L > l/VN + e) < e 



-Ne 



Proof: Recall q = X)i=i <li/N, and define p — J^iLi Pi /N where pi = qi except for i = k, and pk is an arbitrary column 
of the appropriate-sized identity matrix. Then we have 

\\q- <zll 2 - \\P~ 9ll2 < II? -M2 ( fe y triangle inequality) 

i i 

( 1/-ZV) 1 1 q^ — pk 1 1 2 (by definition of p, q and L 2 -norm) 



< (l/JV)Vl 2 + l 2 
= v/2/JV 

This shows that \\q — q\\ 2 is a function of random variables q\, . . . , gjy such changing the fc th random variable q^ for any 
1 < k < N (resulting in \\p — q\\ 2 ) changes the value of the function by at most cj. = y/2/N. Note that g* is not a random 
variable but rather the variable we are trying to estimate. In this case, McDiarmid's inequality (Theorem 19 1 bounds the 
deviation \\q — q\\ 2 from its expectation E \\q — q\\ 2 as: 



Pr(||g — q\\ 2 > E \\q- q\L + e) < exp 



exp 



-2e 2 

-2e 2 
N -2/N 2 



,-JVe 



(7.7) 
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We can bound the expected value using the following inequality: 

2\ V2 



E 



N 



= E 



< E 



JV 



52% 

JV 

E % - N $ 



i=l 



2 / 
2\ 1/2 



1/2 



(by concavity of square root, and Jensens inequality) 



f JV \ 

E E n^-9ii2 

\i=l ) 
t N \ 

^E(£-£) T (£-g) 

U=l / 



1/2 



Multiplying out and using linearity of expectation and properties of $ (namely, that % = 1, E(<fj) = cf and if is constant), 
we get: 

1/2 



E 



JV 



E & - N $ 



/ JV 

^ (E E ( 1 - 2 <^+ Mln < since ^ = 



JV 



JV 



TV 



1/2 



EE(l)-2EE(< ? 1 r 9)+E E ll9l 



vi=l 



i=l 



N-2N\\q\\l + N\\q\\ 2 2 



1/2 



This implies an upper bound on the expected value: 

E\\q~q]\ 2 2 = (l/N 2 )E 



N 



<(1/jV 2 ).jV(1-||^) 
E||^-q1 2 < (l/VN)yJ(l -Ml) 
< (1/VN) 



Using this upper bound in McDiarmids inequality (equation ( |7.7[ )), we get a looser version of the bound that proves the 
proposition: 



Pr(||g-q1| 2 > 1/VN + e) < e 



□ 



We are now ready to prove Lemma 1 8 



PROOF: [Lemma 18 We will treat P\,Pi,\ and P^ Xt i as vectors, and use McDiarmid's inequality to bound the error in 
estimating a distribution over a simplex based on indicator vector samples, using Proposition [20| We know that 



Pr 



Now let r\ = e Ne . This implies 



(\\q-q\\ 2 > Vl/N + e) 

In 77 = -Ne 2 
ln(l/77) = N ^ 
c 



< e 
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Hence, 

Pr q\\ 2 > jl/N + VMV^)/^) <V 

Therefore, with probability at least 1 — 77, 

\\q-q\\ 2 < l/VN+y/]n(Vv)/N 



(7.8) 



Now, in place of q in equation ( |7.8| l, we substitute the stochastic vector Pi to prove the first claim, the vectorized version of 
the stochastic matrix P2,i to prove the second claim, and the vectorized version of the stochastic tensor i-3,2,1 G K™ x ™ x ™ 
obtained by stacking P3 >x ,i matrices over all x, to prove the third claim. The matrices P3 )X ,i are stacked accordingly to 
obtain the estimated tensor £3,2,1. We get the following: 



£i < 1/Vn + y/hi{\ / rf) / N (hence proving the first claim) 
£2.1 < 1/VN + \/\tl{\ / rf) / N (hence proving the second claim) 



3.x,l\i,j [Ps,x,lii,j 



x 1 j 



Pd.,2.1 — Pi 



3,2,1 



Pd.,2,1 — Ps,2,l 

< y/l/N + \/ln(l Jrjj /N (hence proving the third claim) 
Note the following useful inequality from the above proof: 



IJ2 e ^ 2 ^ VW+VHl/ri)/N 



(7.9) 



It remains to prove the fourth claim, regarding e 3,x,i- First we get a bound that depends on n as follows: 

^ £3,x,l = J! l e 3,a;,l | (V Vx, £ 3 ,x.l > 0) 

X X 

< y/H /^e3,x,i 2 (•.■VfeM a , llf^ < ^||x|| 2 ) 



< + J- 



n 1 
In ■ 



N Tj 

We aren't going to use the above bound. Instead, if n is large and TV small, this bound can be improved by removing direct 
dependence on n. Let e(k) be the sum of smallest n — k probabilities of the second observation x^. Let Sk be the set of 
these n — k such observations x, for any k. Therefore, 

e{k) = ]T Pr[z 2 = x] = X v 
Now, first note that we can bound Y^x4S k e 3,x.i as follows: 

x£S k x£S k 

<Vk ea,.,! 2 (•.•VfeIRM|f|| 1 <^||f|| 2 ) 
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By combining with equation 7.9 we get 

E £ 3,x,i < + y/khi(l/v)/N (7.10) 

To bound Ylxes e 3,x.i, we first apply equation ( |7.8| l again. Consider the vector <f of length kn 2 + 1 whose first kn 2 entries 
comprise the elements of Pz.x.i for all x ^ Sk, and whose last entry is the cumulative sum of elements of Pj,. x ,i for all 
x € 5fe. Define q accordingly with P3, x ,i instead of P3, x ,i- Now equation ( |7.8| > directly gives us with probability at least 

1 — rj: 



E / J([- P 3,a;,l]tJ ~ [-P3,x,l]i,j) 2 



E Ea^wk - E E^ 



< y/l/N+y/hi(l/ri)/N 



E E([^,i]« - [Ps 



Since the first term above is positive, we get 

E ECI-^.x.llij - [P3,x,l]ij) 
x£S k i,j 

Now, by definition of Sk'-, 

E^ e3,x,l = E] |-^3,x,l — P3,x,l 
x£S k xeS k 

-EE |[-P3,x,i]y - 

x<ES k i,j 

= E E maX ( ' [-^3^.l]*J ~ [ P 3,a;,l]u 

xes k i,j 



< ^ljN+ y/]n(l/ri)/N 



[P 3 ,x,ih (Wx,p|| 2 < \\x\l,) 



- E E min (°' [-^3,x,i]ij - [Pa,x,i]ij) ('■' Vx, |sf| = [max(0, x) - min(0, £)]) 

xes k i,j 

< E E max (°' [Pa,^h - [ p 3,x,i]«) + E E^aW 

xES k i,j x£S k i,j 

+ E E min (°' [^3,a,l]« ~ /'i < I + E S[ p 3.»,l]« 

= E E max (°> - i p s.x,ih) + <k) 

x£S k i,j 

+ E E min (°' [-^3,x,i]ij - [-P3,a,i]iiJ + e(fc) (by definition of e(fc)) 

E E ( / ''- : v - [*W]«i) + 2e W 

Plugging in equation ( |7.11| i, we get a bound on J2 x es k e 3,x,i'- 

E £ 3 A i < Vl/N + ^ln(l/ V )/N + 2e(k) 

x£S k 



Combining with equation ( |7.10[ ) and noting that k is arbitrary, we get the desired bound: 

E £ 3.M < ™HVkHl/v)/N + y/k/N + ^Hl/ri/N + ^l/N + 2e(k)} 



(7.11) 
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Note that, to get the term ln(3/?/) instead of ln(l/?j) as in the fourth claim, we simply use 77/3 instead of 77. This bound on 
e 3,x.i w iU t> e small if the number of frequently occurring observations is small, even if n itself is large. □ 



The next lemma uses the perturbation bound in Corollary 17 to bound the effect of sampling error on the estimate U, and 
on the conditioning of (U T OR). 

Lemma 21 [Modification of HKZ Lemma 9] Suppose e 2 _\ < e ■ <Jh{Pi,\) for some e < 1/2. Let Eq = e 2 i/((l — 
e)a k (P 2 ^)) 2 . Define U,U € IR mxfc as the matrices of the first k left singular vectors of P 2 .\ 1 P 2 _\ respectively. Let 
9\, ... ,9k be the canonical angles between span(?7) and span(f7). Then: 

1. e < 1 

2. a k (U T P 2A ) > (l-£)a fc (P 2 ,i) 

3. a k (U T P 2A ) > VT^<J k (P 2A ) 

4. a k (U T OR) > ^l^T 0( j k {OR) 

Proof: First some additional definitions and notation. Define J7j_ to be the remaining n — k left singular vectors of P 2 ,i 
corresponding to the lower n — k singular values, and correspondingly U± for P 2 .\- Suppose ITEV T = P 2 ,i is the thin 
SVD of P 2 ,i. Finally, we use the notation ^{A} £ M. q to denote the i th right singular vector of a matrix A 6 M ?)X<? . Recall 
that <7i{A) = \\Avi{A}\\ 2 by definition. 



First claim: So < 1 follows from the assumptions: 



£0 = 



£2.1 



< 



((l-e)a k (P 2 , 1 )f 

e 2 a k (P 2 ,i) 2 
(1 - e) 2 a k (P 2>1 ) 2 

-2 



(1-ef 
< 1 (since e < 1/2) 

Second claim: By Corollary [nj (J k {P 2 .\) > (1 — s)cr k (P 2 ,i). The second claim follows from noting that a k (U T P 2 ,i) 
Third and fourth claims: First consider the fc th singular value of U T U. For any vector lei 1 : 



U T Ux 




U T Uy 




— > min - 





y \\y\\ 2 

a k (U T U) (by definition of smallest singular value) 



cos(#/c) (by Definition 14 1 

^/l-sin 2 (0 fc ) 



'l - a k {Uj_U) 2 (by Lemma[l5 



uju 



(by definition of L 2 matrix norm) 



Therefore, 



U T Ux 



> llxll 



ulu 



(7.12) 
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Note that 



uju 



<e 2A 2 /a k (P 2 ,i) 2 (by Corollary [XT} 



< 



" (l-e) 2 ( r fc (P 2 , 1 ) 
= £o (by definition) 
Hence, by combining the above with equation ( 7.12| i, since < £q < 1: 



(since < e < 1/2) 



U T Ux 



> \\x\Ly/l-e (for all x € K fe ) 



(7.13) 



The remaining claims follow by taking different choices of x in equation ( |7.13| l, and by using the intuition that the smallest 
singular value of a matrix is the smallest possible L 2 norm of a unit-length vector after the matrix has left-multiplied 
that vector, and the particular vector for which this holds is the corresponding right singular vector. For claim 3, let 
x = YV T v k {U T P 2 ,i}. Then by equation ( f7l3) : 



u T uY,v T v k {u T p%i} > zv T v k {u J p 2A } yr^o 



Since P 2A = UT,V T , and ||S^ t j4{S^ t }|| 2 < Y,V T v k {U T P 2<1 } ^ by definition of iJ k {Y,V T }, we have: 

u T p 2A mu t P2.1} > \\zv T is k {zv T }\\ 2 Vi^o 

<J k {U T P 2A ) > <j k {YZV J )VT^T a (by definition of a k {U T P 2A ) , a k {ZV T )) 
a k (U r P 2 ,i) > <?k(P2,iWT^o (■: <r fc (S7 T ) = a k (P 2A )) 

which proves claim 3. 

For claim 4, first recall that OR can be exactly expressed as P 2 ,i(S diag(7?)0 T ) + (equation ( |7.2| i). For brevity, let A = 
(S diag(7?)0 T )+, so that OR = P 2A A. Then, let x = Y,V T Av k {U T OR} in equation (Tl3\ : 



U T UY,V T Av k {U T OR} > T,V r AP k {U T OR} V^eo 
2 2 

Since P 2 1 = UY,V T , and ||SV rT ^i7 fc {EV T ^}|L < Y>V T AV k {U J OR\ by definition of v k {EV T A\, we get: 

11 112 2 

U T P 2 xAM^OR} > \\Y<V T Av k {Y,V T A}L VT^> 
2 11 IM 

U T ORv k {U T OR\ > a k (T l V T A)Vl^ r £o (by equation (fH)) 
2 — 

By definition of a k (U T OR), a k (Y,V T A), we see that 



<r k (u T OR) > a k (j:v T A)Vi^ 

a k (U T OR) > a k {OR)VT^o (V ^V 1 'A) = a k (P 2 ,iA) = <T fc (OiZ)) 



hence proving claim 4. 



□ 



Define the following observable representation using U = U , which constitutes a frae observable representation for the 
HMM as long as (U T OR) is invertible: 

boo = {P2A T U)+P 1 = (U T ORf T R T i m 

B x = (U T P^ xA )(U T P 2A )+ = (U T OR)W x (U T OR)- 1 for x = 1, . .. ,n 
h = U T P 1 



26 



Define the following error measures of estimated parameters with respect to the true observable representation. The error 
vector in Si is projected to W n before applying the vector norm, for convenience in later theorems. 

(U T 0) T (b oo -boo) 

A x = ^OR)- 1 (b x -B x )(U t OR) i= (U t OR)- 1 B x (U t OR)-W x 

X 

Si = jR(U T OR)- 1 (b 1 -bi) = R{U 1 OR)- 1 b l - n 

The next Lemma proves that the estimated parameters boo, B x , bi are close to the true parameters b^, B x , b\ if the sampling 
errors ei, £2.1, £3,2.1 are small: 

Lemma 22 [Modification of HKZ Lemma 10] Assume e 2 \ < ffc(^ > 2,i)/3. Then: 



<5oo <4- 



£2,1 



cr k (P 2 ,i) 2 3a fc (P 2>1 ) 



A, < 



A < 



Pr[x 2 = a;] 



e 2,l 



VS °k(OR) 
8 y/k 

7! ' °k(QR) \a k (P 2tl ) 2 ' 3a fe (P 2ll ) 



^ fe (P 2 ,i) 2 3°k(P2,i) 

£2,1 ^x^3,x,l \ 



ei 



VS <J k (OR) 

Proof: Note that the assumption on £2.1 guarantees (U T OR) to be invertible by Lemma 
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claim 4. 



bound: We first see that 600 can be bounded by 

<5oo 



'00 ^00 



< llo T 

< 



(O^^-ftoo) 

?/(&oo - boo) 



U^br^ - 6^) 
Uiboc - boo) 
boo ^00 



In turn, this leads to the following expression: 



{Pl l U)+P 1 - (P 2 , 1 T U)+P 1 

{p 2 ,iU) + Pi - (p 2 .i t u)+p 1 + (p 2 ,i t ^) + A - (p 2 ,i t c>) + Pi 

((PliU) + - {P 2 ,i J U) + ) Pi + (P 2 ,i T C/) + (Pi - Pi)|| 2 
((PJiU) + - (P 2 ,i T C/) + )Pi + (P 2 ,i T U) + (P 1 -P 1 ) 



((^2 T lf/) + -(^2,l T C/) + ) 



Pi 



(p 2 .l T t/) + (Pi -Pi) 

2 
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The last step above obtains from the consistency of the L 2 matrix norm with L\ vector norm, and from the definition of L 2 
matrix norm (spectral norm) as ||vl|| 2 = max ^pjp ■ Now, recall that U has orthonormal columns, and hence multiplying a 
matrix with U cannot increase its spectral norm. Hence, 

P 2 T 1 U-P 2 ,i T U = (Pj, 1 -P2,i T )U < Pl 1 ~P 2 ,i T 



£2,1- 



So, we can use Lemma 
bound on the differenc 
trix are the reciprocals 

^k(Pj A U) = a k {P 2A ) 



12 



to bound the L2-distance between pseudoinverses of PjiU and P 2 ,i T U using £2.1 as an upper 
bound on the difference~between the matrices themselves. Also recall that singular values of the pseudoinverse of a ma- 
trix are the reciprocals of the matrix singular values. Substituting this in the above expression, along with the facts that 



Pi 



= 1 and 



< 



(Pi -Pi) 

1 + V5 



ei, gives us: 



£20 £1 

2 mm(a k (P 2A ),a k (P 2A T U)Y *k(Pi,l T U) 



Now, to simplify the last expression further, consider Lemma 21 in the above context. Here, £2,1 < 0fc(P 2i i)/3 and hence 



e = 1/3. Therefore a k (P 2 ,i) = <r k (U T P 2A ) > (2/3)a k (P 2A ) and a k (U T P 2tl ) > V±^oVk(P 2 ,i)- Hence 

min (o k {P 2 ,i),a k {P 2 s T U)f = cr k {P 2A f • min(2/3, x/T^) 2 
The latter term is larger since 



£o 



£2,1 



{{l-e)a k {P 2A )Y 



< 



v k (P 2 ,i) /Q 



~ 4a fe (P 2>1 )79 
-1/4 

=> VT^s^ > V3/2 > 2/3 

Therefore min (cr k (P 2 .i), < J k(P 2 ,i T U) S J > o k (P 2 .i) 2 (2/3) 2 . Plugging this into the expression above along with the fact 
that <J k (U T P 2 ,i) > (v / 3/2)crfe(P 2i i), we prove the required result for Soc,: 

1 + V5 9e 2 ,i 2ex 



< 4 



2 4a fe (P 2 ,i) 2 V3<r k (P 2 ,i) 
( ^ + 61 ) 

\<Tk(P 2 ,l) 2 °h{Pi,l)) 



A X ,A bounds: . We first bound each term A x by \/k 



B T — B a 



/a k ((U T OR)): 



A x = {U T OR)- 1 [B x - B x j (U T OR) 

< (U T ORy 1 (B x - B X )U T \\OR\\-y (by norm consistency) 

< Vk (U T ORy 1 (B x - B X )U T WOR^ (by L 1 vs. L 2 norm inequality) 

lOUx ||P||i (by norm consistency) 



< y/k 



(U T ORY 



< y/k (U T OR)- 1 



2 


B x 


— B x 


to 




2 


B x 


— B x 


2 


( 



2 



Vk 



B x — B a 



\o\\ x ,\\R\\ x <i) 

/a k (U T OR) (■: a^fFOB)- 1 = l/a min (U T OR)) 
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The term 



B x - B x 
B x — B 



in the numerator can be bounded by 



< 



(U T P 3 , x ,i)(U T P 2A ) + - (U T P 3 ^ 1 )(U T P 2A ) + 
(U T P 3 , xA ) ((U T P 2A )+ - (f> T P 2 ,i) + ) || 2 + ||^ T (•Pa.x.i - P 3 ,x,i) (U T P 2A )+ 



< \\P 



1 + V5 



£ 2,1 



3,x,l\\ 2 



< Pr[x 2 = x] 



min (a k {P 2 ,i),a k (U T P 2 ,i) 
1 + \/5 e 2 ,i 



e 3,x,l 

<J k (U J P 2A ) 



iin (o*(ft,i), ^(t/ T P 2 ,i)) ^(^2,1) 



where the second inequality is from Lemma 12 and the last one uses the fact that 



liWlla < II^WII* = . /E[ p 3-,i]L- ^ X ^ * ' ' = Pr ^ = 



Applying Lemma 21 as in the 6^ bound above, gives us the required result on A x . Summing both sides over x results in 
the required bound on A. 

5\ bound: For Si, we invoke Condition |4]to use the fact that ||-R||j < 1. Specifically, 



R(U T OR)- 1 U T (P 1 - P t ) 



< \\R\ 



(U T OR)- 1 U T (P 1 - Pi) 



(norm consistency) 



^VkWR^ (U T OR)- 1 U T (P 1 ~ Pi) Cll^lli < Vn|N| a foraiiya:GM n ) 



KVkWR]^ {U T OR)- i: U J (Pi -A) 



(norm consistency) 



< Vk \\R\\ 1 (U J ORy l ■ ei (defn. of e x , U T has orthogonal columns) 



v k (U T OR) 



\R\li < 1, defn. of L 2 -noim) 



The desired bound on Si is obtained by using Lemma 21 With e, Eq as described in the above proof for 6^, we have that 
Uk{U T OR) > (y/3/2)<Jk(U T OR). The required bound follows by plugging this inequality into the above upper bound for 
Si. " ' □ 



7.4 Proof of Theorem |2] 



The following Lemmas [23] and 24 together with Lemmas 18|21|22 above, constitute the proof of Theorem [2] on joint 



probability accuracy. We state the results based on appropriate modifications of HKZ, and provide complete proofs. We 
also describe how the proofs generalize to the case of handling continuous observations using Kernel Density Estimation 
(KDE). First, define the following as in HKZ 

e(i) = min i ^Pr[x 2 = j] : S C {1 . . . n}, \S\ = li- 
lies 

and let 

no(e) = min{« : e(i) < e} 

The term no(e), which occurs in the theorem statement, can be interpreted as the minimum number of discrete observations 
that accounts for 1 — e of total marginal observation probability mass. Since this can be much lower than (and independent 
of) n in many applications, the analysis of HKZ is able to use n instead of n in the sample complexity bound. This is 
useful in domains with large n, and our relaxation of HKZ preserves this advantageous property. 

The following lemma quantifies how estimation errors accumulate while computing the joint probability of a length t 
sequence, due to errors in B x and b. 



29 



Lemma 23 [Modification ofHKZ Lemma 11] Assume U T OR is invertible. For any time t: 

]T jR^OR)- 1 (b x Ji - B Xt -M) ^ < (1 + A)*Ji + (1 + A)* - 1 



Proof: Proof by induction. The base case for t — 0, i.e. that 



i?(C/ T Oi?)- 1 (&i-6i) 



< Si is true by definition of 



Si. For the rest, define unnormalized states b t — &t(£i:t-i) = Sx t -i:i&i an d &t = &t(^i:t-i) = £a: t _i:i&i. For some 
particular t > 1, assume the inductive hypothesis as follows 

^ |i?(C/ T 0i?)- 1 (fe t - &t) J < (1 + A)*5i + (1 + A)' - 1 

Hit 

The sum over xi : t in the LHS can be decomposed as: 
2|ii(tf r Oi2)- 1 (& t -& t )| ] 

= E E I^^OiZ)- 1 ((Sx t - B xt )bt + (B xt - B xt )(bt -bt) + B xt (b t -bt)) 

X Xl;t-1 

Using triangle inequality, the above sum is bounded by 

£ ]T \\r(U t OR)- 1 [B xt - B xt ) (U T O)\\ i \\R0 T OR)- 1 b t \\ i 

+ E E \\R0 T OR)- 1 ( K B Xt -B Xt )(U T O)\\ i \\R(U T OR)- 1 (b t -b t 

X t Xl-.t-l 

+ E E \\R(U T OR)- 1 B t (U T OR)(U T OR)- 1 (b t -b t ) 



Each of the above double sums is bounded separately. For the first, we note that 



R(U T OR)- 1 b t = Pr[xi :t _i], which 



sums to 1 over Xi-.t-i- The remainder of the double sum is bounded by A, by definition. For the second double sum, the 



inner sum over 



Ri^ORj-^bt-bt) 



is bounded using the inductive hypothesis. The outer sum scales this bound by 



A, by definition. Hence the second double sum is bounded by A((l + A)*^ 1 ^! + (1 + A)* -1 — 1). Finally, we deal with 
the third double sum as follows. We first replace (U T OR)' 1 B t (U T OR) by W Xt , and note that R ■ W Xt = A Xt R. Since 
A Xt is entry-wise nonnegative by definition, || A^v^ < 1^^4 Xt |tT|, where \v\ denotes element-wise absolute value. Also 
note that \v\ = t^T\v\ = = \\v\\ v Using this result with v = R{U T OR)~ l (b t ~b t ) in the third double 

sum above, the inductive hypothesis bounds the double sum by (1 + A)* -1 ^ + (1 + A)* -1 — 1. Combining these three 
bounds gives us the required result: 



J2\\R(u T OR)- 1 (b t -b t ) 



< A + A((l + A)*- 1 ^ + (1 + A)*- 1 - 1) + (1 + A)'- 1 ^ + (1 + A)*- 1 
= A + (1 + A)((l + A)'- 1 ^ + (1 + A)*- 1 - 1) 
= A + (1 + A)*5i + (1 + A)*-1- A 
= (1 + A)*Ji + (1 + A)* - 1 
thus completing the induction. 

The following lemma bounds the effect of errors in the normalizer b^. 

Lemma 24 [Modification of HKZ Lemma 12] Assume e 2 ,i < Cfe(P2,i)/3. Then for any t, 



□ 



Pr[a;i :t ] - Pr[x 1:t ] < (1 + <5 co )(l + Si)(l + A)* - 1 
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Proof: First note that the upper bound on 62.1 along with Lemma \2 1 1 ensure that a k (U T OR) > and so (U T OR) is 
invertible. The LHS above can be decomposed into three sums that aredealt with separately: 

]T |Pr[x 1:t ] - Pr[x 1:t ]\ = J2 \blB Xt:1 h - b^B^ J, 

Xu.t Xut 

Xl: t 

< J2 \(boo-b oo ) J {U J 0R){U J 0R)- 1 B Xt ~b 
+ J2 \(boo~b 00 ) T {U T OR){U T OR)-\B Xt l b 1 - B Xt . A bi) 

Xut 

+ E ^li^OR^ORr'iB^h - B xt:1 h) 

xut 

The first sum can be bounded as follows, using Holders inequality and bounds from Lemma [22] 

Y I (600 - b 00 ) T (U T OR)(U T OR)- 1 B Xt l b 1 \ < Y ||(C/ T 0) T (6 00 - &«,) R(U T OR)~ 1 B Xt .J>i 
* — ' l 1 A — ' ll 00 

Xl;t Xl-t 

< y$oo Px t:1 7fiii 

Xut 

= y 800 Pr[a;i:t] 

Xu.t 

= Soo 

The second sum can be bounded also using Holders, as well as the bound in Lemma [23] 

y I (boo ~b 00 ) T {U T OR){U T OR)-\B Xt . 1 b 1 - B Xul bx) 

Xl:t 

< (U T 0) T (b oc -b oc ) R{U T OR)-\B Xtl b 1 -B Xt . 1 b 1 ) 

OO 

<5 00 ((l + A) t J 1 + (l + A) t -l) 
The third sum again uses Lemma [23] 

y \bJ a (U T OR)(U T OR)- 1 (B Xt . 1 b 1 - B Xt ~b x )\ = y \l Y R{U T OR)- 1 {B Xt . A b l - B^Ji) 

Xl.t Xut 

< \\R(U T OR)- 1 (B Xt:1 b 1 - 5^,61)^ 

< (1 + A)*$i + (1 + A)* - 1 

Adding these three sums gives us: 

y |Pr[x 1:t ] - Pi[x la ] \<S QO + ^((1 + A)**! + (1 + A)* - 1) + (1 + A)% + (1 + A) ( - 1 

xut 

< 5^ + (1 + <5 co )((l + A)**i + (1 + A)* - 1) 
which is the required bound. □ 
PROOF: (Theorem[2]). Assume N and e as in the theorem statement: 



e = a k (OR)a k {P 2 , 1 )e/(UVk) 
N>C t i ( 



k ■ n (e) 



e 2 \o k {ORYa k {P 2tl Y ( r fc (Oi?) 2 a fc (P 2)1 ) 2 



log(l/77) 
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First note that 



^2\Pr[x 1:t ]-Pi[x 1:t ] 



< 2 



since it is the L\ difference between two stochastic vectors. Therefore, the theorem is vacuous for e > 2. Hence we can 
assume 

e < 1 

in the proof and let the constant C absorb the factor 4 difference due to the 1 /e 2 term in the expression for N. 
The proof has three steps. We first list these steps then prove them below. 



First step: for a suitable constant C, the following sampling error bounds follow from Lemma 18 



ei < min (.05 • (3/8) • a k {P 2 ,i) ■ e, .05 • (^3/2) • a k (OR) ■ (l/Vk) • e) 
e 2 ,i < min (.05 • (1/8) • a k (P 2 ,i) 2 • (e/5), .01 • (V5/8) ■ a k (OR) ■ a k (P 2 .i) 2 • (l/(tV*)) ■ e) 
£ 3,-,i < 0-39 • (3V3/8) • a k {OR) ■ a k (P 2A ) ■ (l/(Wk)) • e 



(7.14a) 
(7.14b) 
(7.14c) 



Second step: Lemma 22 together with equations ( |7.14| i imply: 



(Soo < -05e 
Si < .05e 
A < OAe/t 



Third step: By Lemma 24 equations ( |7.15[ ) and the inequality 

(1 + (a/t)Y < 1 + 2a for a < 1/2 

we get the theorem statement. 

Proof of first step: Note that for any value of matrix P 2 ,i, we can upper-bound a k (P 2y i) by 1: 

1k(P 2 ,l) < <Tl{P 2 ,l) 



max ||iVix|| 



IMI,=i 



n / n 



1/2 



,iE E E^]^ 



j=l \i=l 



< max 



2 ,l\ijX, 



(by norm inequality) 



^ EE l[ p 2,i]yl (N < 1 since ||z|| 2 = 1) 

j=l i=l 
n n 

= EE^« ( by non - ne g ativit y of p 2 , i ) 

= 1 (by definition) 



(7.15a) 
(7.15b) 
(7.15c) 



(7.16) 
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Similarly, for any column-stochastic observation probability matrix O we can bound (OR) by \fk. First see that a\ (O) < 



o"i(0) = max ||0a;|| 2 

1 1 x I L — 1 



1/2 

\2 



max yzyziOijXiY 

x L=l I — ~ 

/ \ 1/: 

/ m n \ 

! EE ?? (N| 2 = l=i-N<l) 



m n 



< max 

1/2 

- I y^(y^Qtf) 2 1 (by triangle inequality) 

/ m \ V2 

< hTl 2 (by definition of O) 



'to: 



Now the bound on a k (OR) follows from Conditional. e. a k (OR) < \Jk/r 
a k (OR) = min \\ORx\\ 2 

l=c|| 2 = 1 

< ||0|| 2 • min ||i?a;|j 2 (by norm consistency) 



\ x IU— i 



< ym min 

I a; I L — 1 



^2 y^X^io^) 2 ('•' 11^112 = for an y matrix A) 



Assume the c*' 1 column of i? obeys Condition p\ for some 1 < c < k. Also assume a: = e c , the c th column of the k x k 
identity matrix, which obeys the constraint ||x|"= 1. Then every component of the inner sum is zero except when j = c, 
and the min expression can only get larger: 



\ i=l 

= Vm\\R[;c]\\ 2 
= \fk 

hence proving that a k (OR) < \[k. 

Now we begin the proof with the ei case. Choose a C that satisfies all previous bounds and also obeys (vC/4)-0.05-(3/8) > 
1. 

ei < y/l/NiVHS/v) + 1) (by Lemma[l8]l (7.17) 
< y/l/N(2y/ln(3/r])) (since a/M 3 /^) > \/hi3 > 1) (7.18) 
Now, plugging in the assumed value of 7Y: 



£i < 2e(a k (P 2 , 1 ) 2 a k (OR)) \n(3/ V ) (? ig) 



t^Ck(l + n (e)a k (P 2A ) 2 ) V ln (l/*?) 
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Any substitutions that increase the right hand side of the above inequality preserve the inequality. We now drop the additive 
1 in the denominator, replace y/\a{Z/rj) / ln(l/^) by 2 since it is at most Vln 3, and drop the factors t, J rio(e) from the 
denominator. 



ei < 



4a k (OR)a k (P 2A )e 



^Cka k (P 2>1 ) 



(7.20) 



1 



< 4 
1 



a k {OR)/^k\ ■ [a k (P 2 ,i)]-e 
< -^=min(4-<7 fe (P 2 ,i) ■e,4-a k {OR) ■ • e) (.both [a fc (Oi?)/%/fc] and [a k (P 2 ,i)] are < 1) 



^min (0.05 • 3/8 • <r k (P 2 ,i) ■ e,0.05 ■ y/3/2 ■ <J k (OR) ■ 1/Vk • e) 



(forC = • 0.05 • 3/8) 



= min (o.05 -3/8 ■ a k (P 2 ,i) ■ e,0.05 • V3/2 ■ cr k (OR) ■ l/Vk-e\ (v C > 1) 
Hence proving the required bound for t\. 



Next we prove the e 2 ,i case. Choose a C that satisfies all previous bounds also obeys (yC/4) ■ 0.01 • (V3/8) > 1. Note 
that, since the bound on e 2i i in Lemma 18 is the same as for €%, we can start with the analogue of equation ( |7.19| i: 



£2,1 < 



2e{a k {P 2 ,ifa k {OR)) H*/v) 



tJCk{\ + n (e)a fe (P 2il ) 2 ) V ln (V?7) 



We now drop the additive no(e)a k (P 2t i) 2 in the denominator, again replace -\/ln(3 /rf)/ ln(l /?y) by 2 since it is at most 
Vln 3, and drop the multiplicative factor t from the denominator. 







< 




















< 


1 


— = mm 




y/C 


< 


1 . / 


— mm 




c \ 



(T k (P 2 ,lY 



cr k 



(OR)/Vk 



<D 



(for C = (VC/4) • 0.01 • (V3/8)) 
< min (o.05 • 1/8 • cr fe (P 24 ) 2 ■ e,0.01 • V3/8 • a k {OR) ■ cr k (P 2A f ■ l/Vk-e) (since C > 1) 



2-0.39-(3y / 3/8)v / C 



hence proving the bound on e 2 



Finally for J2 X e 3,x,i, assume C such that — 
the bound from Lemma [T8l 



> 1 in addition to previous requirements on C. we first restate 



£ 3,-i < mm (-/JIN (y/fo3/fj + l) + 2e(j)) + y/l/N (y/hiS/n + l) 

X 

< ^n (e)/N (Vln3/?7 + l) + 2e(no(e)) + ^l^V (^3^? + 1 

< v^lTiV (v/ln3/?7 + l) (no(e) + 1) + 2e (since e(n„(e)) < e) 
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The first two terms are exactly as before, so we perform the same steps as in equations ( |7.17| i-( |T20] i except we do not drop 

ty/no(e), to get: 



e 3,*,i < ln , . I («o(e) + 1) + 2e 

< 4 ^(°^( P ^) £ . ( 2 . Me)) + 2a k (OR)a k (P 2il )e/4tVk 

ty/CKn {e) 

(since 1 + n (e) < 2 ■ n (e), and plugging in e) 

< <j k (OR) ■ cr k (P 2 ,i) ■ tVk ■ e ■ (8/VC+1/2) 



< ^0.39 • (3V3/8) • a k (OR) ■ <7 fe (P 2jl ) • tVfc • e (for C = ^gg^l^g ) 

< 0.39 • (3^3/8) • a k (OR) ■ «7 fc (P 2 ,i) • ■ e (since C" > 1 by assumption) 



Hence proving the required bound for J2 X e 3,x,i- 

Proof of second step: Substituting from equation ( |7.14| > into Si in Lemma 
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2 

< .05e 



Substituting from equation ( |7.14| i into Soo in Lemma 22 



See <4 



£2,1 ei 



^■(-P 2 ,i) 3tr fc (P 2 ,i) y 
< / x2 min( .05 • (1/8) • cr fc (P 2 ,i) 2 • (e/5), .01 • (V3/8) • a k (OR) • <7 fc (P 2)1 ) 2 • (l/(tVk)) ■ 

CTfc(P 2 .l) V 

05 • (3/8) ■ (7 fe (P 2 i) ■ e, .05 ■ (73/2) • a fc (OP) ■ (1/%/fc) 



o-fc(P 2> ij 
4 



mm 



3a fe (P 2 ,i) 

< min( .05e, .04 • (V3/8) • a k (OR) ■ (l/(tVk)) 



+ minf .05 • (1/2) • e, .05 ■ (2/V3) ■ ^TT^T ' (V^) ' ^ 
V cr k (F 2 ,l) ) 



< .05e(.01 + .5) 

< .05e 
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Substituting from equation ( |7.14| i into A in Lemma 22 

\fk ( £2,1 ^x^3,x,l 



A < 



< 



^3 <T k {OR) \a k {P 2A f 3o- fc (P 2 ,i) 
&Vk ( 1 



V3(T fc (Oi?) V ff fc( P 2,l) 

1 



(.05 • (1/8) • a fe (P 2 ,i) 2 • (e/5), .01 • (v^/8) • a k (OR) ■ a k {P 2A ) 2 ■ 



3f7 fe (P 2a ) 
min .05 • (e/5) 



0.39 • (3V3/8) • <J k (OR) ■ a k (P 2A ) • 



^ oi.5 

V3a fc (Ofl)' * 



0.39 



< .01 • - +0.39 • - 

t t 



< OAe/t 
Proof of third step: By Lemma 
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5^ |Pr[a:irt] - Pr[si !t ] < (1 + *«,)(! + <5i)(l + A)* - 1 



< (1 + .05e)(l + .05e)(l + 0.4e/t)* - 1 (by equations ( fTB) ) 

< (1 + .05e)(l + .05e)(l + 0.8e) - 1 (by equation ( p77T6] >, since 0.4e < 1/2) 
= 1 + .05e + .05e + .05 2 e + 0.8e + .04e 2 + .04e 2 + (.05) 2 • .08e 3 - 1 

= .0002e 3 + .0825e 2 + 0.9e 

< (.0002 + .0825 + 0.9)e (since e < 1 by assumption) 
= 0.9827e 

< e 

This completes the proof of Theorem|2] 

7.5 Proof of Theorem [2] for Continuous Observations 



□ 



For continuous observations, we use Kernel Density Estimation (KDE) [ 19] to model the observation probability density 
function (PDF). We use a fraction of the training data points as kernel centers, placing one multivariate Gaussian kernel 
at each point. 1 The KDE estimator of the observation PDF is a convex combination of these kernels; since each kernel 
integrates to 1, this estimator also integrates to 1. KDE theory [19| tells us that as the number of kernel centers and the 
number of samples go to infinity and the kernel bandwidth goes to zero (at appropriate rates), the KDE estimator converges 
to the observation PDF in L\ norm. The kernel density estimator is completely determined by the normalized vector of 
kernel weights; therefore, if we can estimate this vector accurately, our estimate will converge to the observation PDF as 
well. 

Hence our goal is to predict the correct expected value of this normalized kernel vector given all past observations (or more 
precisely, given the appropriate sequence of past observations, or the appropriate indicative events/features). In the context 
of Theorem [2] joint probability estimates for i-length observation sequences are effectively the expectation of entries in a 
t-dimensional tensor formed by the outer product of t indicator vectors. When we move to KDE, we instead estimate the 
expected outer product of t stochastic vectors, namely, the normalized kernel weights at each time step. As long as the sum 
of errors in estimating entries of this table goes to zero for any fixed t as the number of samples increases, our estimated 
observation PDFs will have bounded error. 



The only differences in the proof are as follows. In Lemma 18 we observe to be stochastic vectors instead of indicator 



vectors; their expectation is still the true value of the quantity we are trying to predict, pi are also stochastic vectors in that 
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proof. In the proof of Proposition 
always equal to 1, and the same holds for pjpi. Also 
things, the above proof goes through as it is 



p k is an arbitrary stochastic vector. Also 



< 



= 1 now instead of being 



\Pi — Pi I 



< 



\\Pi — Pi Hi = 1 (by triangle inequality). Besides these 



1 We use a general elliptical covariance matrix, chosen by SVD: that is, we use a spherical covariance after projecting onto the singular 
vectors and scaling by the square roots of the singular values. 
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Note that in the continuous observation case, there are continuously many observable operators W x that can be computed. 
We compute one base operator for each kernel center, and use convex combinations of these base operators to compute 
observable operators as needed. 



8 Appendix II: An Example of Learning with Ambiguous Observations 

When stacking observations, the modified, larger P 2t i € R nxn still has rank at most k since it can be written in the form 
^2,1 = GTH for some matrices G,H T e M" xm . For example, if n = 2 for an HMM with ambiguous observations, 
and we believe stacking 2 observations per timestep will yield a sufficiently informative observation, the new observation 
space will consist of all n = n 2 = 4 possible tuples of single observations and P 2 ,i G K" xn , with each observation i 
corresponding to a tuple < i\, i 2 > of the original observations. Specifically, 

-P2,iCM) = Pr(x4 = j 2l x 3 = j u x 2 = ii,x\ = h) 

= ^2 Pr(x 4 = j 2 , x 3 = j lt x 2 = i2, x 1 =i 1 ,h 4 = d, h 3 = c,h 2 = b, h x = a) 

a.b,c.d 

^ ^ Gj 2 d r ^ 1 dcGj lC T c f ) Oi 2 f ) T} )a Oi ia 7r a 

a.b,c.d 

= °3> cTch \- diag(7r)0 T ] 6 , l where 0^ c = ^ O hd T dc O nc 

b.c d 

^T\i = OT diag(7r)0 T 

Similarly, we can show that P 3 , x ,i = GTH T for some matrices G,H T e E" xm . The exact formulae will differ for 
different choices of past and future observable statistics. 



9 Appendix III: Synthetic Example RR-HMM Parameters 



Example 1 



0.3894 0.2371 0.3735 
0.2371 0.4985 0.2644 
0.3735 0.2644 0.3621 



O 



0.6000 0.2000 0.2000 
0.2000 0.6000 0.2000 
0.2000 0.2000 0.6000 



(9.1) 



Example 2 



0.6736 0.0051 0.1639 
0.0330 0.8203 0.2577 
0.2935 0.1746 0.5784 



O 



1 .5 
1 .5 



Example 3 



0.7829 0.1036 0.0399 0.0736 

0.1036 0.4237 0.4262 0.0465 

0.0399 0.4262 0.4380 0.0959 

0.0736 0.0465 0.0959 0.7840 



O = 



10 10 
10 1 
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